JOURNAL OF PLANKTON RESEARCH

j

VOLUME

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation W. C. GENTLEMAN* AND A. B. NEUHEIMER DEPARTMENT OF ENGINEERING MATHEMATICS AND INTERNETWORKING, DALHOUSIE UNIVERSITY, CANADA B3J

1340 BARRINGTON ST, HALIFAX, NOVA SCOTIA,

1Y9

*CORRESPONDING AUTHOR: [email protected] Received April 30, 2008; accepted in principle July 18, 2008; accepted for publication July 23, 2008; published online July 25, 2008 Corresponding editor: Roger Harris

Modellers have long been aware that the mathematical form of zooplankton mortality, or closure, significantly affects the dynamics of planktonic ecosystem models. Another important formulation is the functional response, i.e. how ingestion rates change with prey density. Here we explain why different grazing responses can have profoundly differing influences on modelled dynamics, and how common practices may limit models due to misguided characterization of feeding behaviours. Use of different ingestion functions in a Nutrient – Phytoplankton – Zooplankton (NPZ) model results in oscillating versus constant densities. Contrary to the conclusions of previous studies, it is shown that these results are not due to zooplankton satiation versus non-satiation. Analysis of a predator-prey model is used to derive the necessary condition for ecological stability, which is related to food-limited clearance rates. Sensitivity studies demonstrate that zooplankton clearance rates have a strong influence on the dynamics of more complex models. Moreover, it is shown that acclimation time lags can dramatically alter results from those where zooplankton instantly adapt to changing prey densities due to the corollary effect on clearance rates. These results are discussed in terms of practical advice to modellers who face uncertainty in choosing expressions for the functional response.

I N T RO D U C T I O N Nutrient – Phytoplankton – Zooplankton (NPZ) models and their variants (e.g. NNPPZZ or NPZD, where D is Detritus) have been used for decades to study the dynamics and productivity of planktonic ecosystems (e.g. Steele, 1974; Wroblewski, 1977; Evans and Parslow, 1985; Frost, 1987; Frost and Franzen, 1992; Franks and Chen, 1996, 2001; Edwards et al., 2000a, b; Leising et al, 2003), investigate biogeochemical cycling (e.g. Wroblewski et al., 1988; Fasham et al., 1990; Doney et al., 1996; Denman and Pen˜a, 1999, 2002; Hood et al., 2003) and specify prey fields for models of copepods and fish (e.g. Carlotti and Wolf, 1998; Batchelder et al.,

2003; Aydin et al., 2005; Hinckley et al., 2007). These models are framed as sets of ordinary differential equations describing the rates of change of biomass and density, due to processes such as photosynthesis, respiration, grazing, egestion and re-mineralization. They are often coupled to physics through light and temperature effects on physiological rates, as well as circulation effects on distributions. Their ability to simulate wide ranging environments is due to the wide ranging dynamics arising from different model formulations. It has long been recognized that the mathematical form of zooplankton mortality, or “closure”, can

doi:10.1093/plankt/fbn078, available online at www.plankt.oxfordjournals.org # The Author 2008. Published by Oxford University Press. All rights reserved. For permissions, please email: [email protected]

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

dramatically alter model results. The simple choice of linear versus quadratic closure can cause the model to produce oscillations (e.g. predator-prey cycles) versus constant densities (Steele and Henderson, 1981, 1992; 1995; Murray and Parslow, 1999; Edwards and Yool, 2000, Neubert et al., 2004). Another important formulation is the functional response, which characterizes how zooplankton ingestion rates change with phytoplankton density. Numerous studies suggest that the functional response has a large influence on modelled dynamics (e.g. Holling, 1965; Steele, 1974; Oaten and Murdoch, 1975a, b; Evans, 1977, Myerscough et al., 1996; Gentleman et al., 2003; Fussman and Blasius, 2005), including the classic paper by Franks et al. (Franks et al., 1986a; hereafter referred to as FWF), which demonstrated that the use of two different functional responses resulted in oscillating versus constant densities in a manner similar to different closure schemes. The two responses used in FWF characterized contrasting grazing behaviours. The first was a common satiating function, wherein ingestion rates are effectively independent of prey density once a critical value is reached. The second was a non-satiating function representative of herbivore acclimation to high phytoplankton densities (Mayzaud and Poulet, 1978). The authors’ conclusion that satiating responses lead to oscillations whereas non-satiating responses damp oscillations has been re-iterated in modern reviews (e.g. Davidson, 1996; Franks, 2002). It has also motivated the use of the non-satiating Mayzaud – Poulet response by more than 25 other model applications (e.g. Franks et al., 1986b, Wroblewski and Richman, 1987; Wroblewski, 1989; Doney et al., 1996; McClain et al., 1996; Keen et al., 1997; Leonard et al., 1999; Botte and Kay, 2000; Friedrichs and Hofmann, 2001; Christian et al., 2002; Le Fouest et al., 2005; Salihoglu and Hofmann, 2007), many of which are cited dozens of times, including the second most-cited ecosystem modelling paper (Arhonditsis et al., 2006). FWF also investigated the influence of temporal lags in the functional response due to zooplankton having to acclimate to changes in prey density. Specifically, FWF, altered a traditional NPZ model wherein changes in phytoplankton density instantaneously affect ingestion rates, to one where there is a delay between changes in food and changes in grazing. When incorporated in other processes, such as mortality or re-mineralization, such time lags have been shown to introduce oscillations into otherwise stable dynamics (e.g. May, 1973; Steele and Henderson, 1981; Ruan, 2001). However, FWF found only a minimal effect of acclimation lags, and we know of no modern application that accounts for them.

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

Here, the NPZ model of FWF is used as a template to demonstrate the pivotal role of the functional response on planktonic ecosystems. We show that the choice of formulation, including seemingly subtle differences among satiating responses and time lags, can have a profound effect on model results, and that common beliefs about how feeding behaviours influence dynamics may be misguided. We also present a mathematical condition governing ecological stability, which is interpreted biologically and provides a basis for understanding the dynamical influence of satiation, food-limitation and acclimation. These results are discussed in terms of practical advice to modellers who must face uncertainty in choosing formulations.

METHODS The original NPZ model of FWF The classic NPZ model of FWF, which is used as the basic framework for our investigations into the role of the functional response, is presented here with some notational changes for clarity: dP N ¼ m P  IZ  mP P dt kN þ N max

ð1Þ

dZ ¼ aIZ  mZ Z dt

ð2Þ

dN N ¼ m P þ ð1  aÞIZ þ mP P dt kN þ N max þ mZ Z

ð3Þ

The rate of change of phytoplankton biomass [P, equation (1)] is due to the balance between growth, grazing and cell death. Growth is nutrient-limited according to Michaelis–Menten kinetics (Michaelis and Menten, 1913) also known as Monod kinetics (Monod, 1942) with a half-saturation value kN and a maximum specific growth rate mmax. Grazing depends on I, the zooplankton specific ingestion rate, which itself depends on phytoplankton density according to the functional response (discussed in later sections). Cell death is characterized by the mortality rate mP , which is the fraction of the phytoplankton community that dies per day. The rate of change of zooplankton biomass [Z, equation (2)] is due to the balance between growth and mortality. Growth is a fraction, a, of ingestion, where a is the assimilation efficiency. Mortality is represented by a linear closure term, mZZ, where mZ is the fraction of the

1216

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Table I: Parameters for FWF model mmax kN mP a mZ P0 Z0 N0 a [kP] Imax b c

2 day21 1 mmol N m23 0.1 day21 0.7 0.2 day21 0.2 mmol N m23 0.2 mmol N m23 1.6 mmol N m23 1 [mmol N m23]21 [.7 mmol N m23] 1.5 day21 1 (mmol N m23)21 1.5 (mmol N m23)21 day21

Phytoplankton maximum specific growth rate Nutrient uptake half-saturation constant Phytoplankton mortality rate Assimilation efficiency Zooplankton mortality rate Phytoplankton density initial condition Zooplankton density initial condition Nutrient density initial condition Ivlev decay constant [Corresponding half-saturation constant] Ivlev maximum specific grazing rate Mayzaud –Poulet decay constant Mayzaud –Poulet scaling constant

zooplankton lost per day. In FWF, mortality is attributed to predation, but such a formulation could equally include natural death and/or physical losses (e.g. Frost and Franzen, 1992). Changes in nutrient density [N, equation (3)] are controlled by the characterization of this model as a conservative system such that total mass (P+N+Z) remains constant, and dP/dt+dZ/dt +dN/ dt=0. Nutrient density is thus decreased by uptake from the phytoplankton, and increased by re-mineralization of dead phytoplankton and zooplankton. The parameter values used in the original model (Table I) were based on typical data for copepod-diatom ecosystems. In the two decades since this model was introduced, it has been used in dozens of other studies, including applications for macro-, meso- and microzooplankton. As such, the functional response parameters have ranged in value by more than an order of magnitude. Our simulations use the original FWF model and modifications to it that include variations in both the parameters and the functional forms of I, as well as other terms (e.g. closure). These modifications are discussed in the relevant sections below.

differential equations derived from a Taylor series expansion (Nayfeh and Mook, 1979). This set of equations can be written in matrix form dDX/dt = JDX, where DX = fDP, DZ,. . .g and J is the Jacobian matrix (also known as the “community matrix” in ecology, May, 1973), whose elements are partial derivatives of the original model equations (e.g. @(dP/dt)/@P, @(dZ/dt)/@P, etc.) evaluated at steady state. The solution is a sum of exponential terms, where the exponents equal the real parts of the Jacobian’s eigenvalues (l ), which are the roots of the characteristic equation (Strang, 1988). Thus, when the real parts of all eigenvalues are negative, deviations will decay and the steady state is stable. However, when the real part of even one eigenvalue is positive, deviations will grow and the steady state is unstable. If none are positive, but at least one has a real part equal to zero, the system stability must be determined by higher-order terms. For a predator-prey system described by the rates of change of phytoplankton and zooplankton (PZ), the corresponding Jacobian matrix is  

J¼

Stability analyses of a PZ ( predator-prey) model In addition to performing time series simulations with NPZ models, our study employs stability analyses to investigate the influence of the functional response. Such analyses are centred on the systems’ steady states, which are those variable values for which the source and loss rates are exactly balanced (i.e. values of P and Z for which dP/dt = dZ/dt = 0 etc.). Characterization of steady states is based on how the system behaves when variables deviate slightly from equilibrium values: they are stable when variables return to their steady state, or unstable when deviations grow. This qualitative distinction is made by considering the terms governing the first-order time variation of the deviations from steady state (DP(t), DZ(t), etc.) as described by a set of linear

@ dP @P dt   @ dZ @P dt

 !

@ dP @Z dt   @ dZ @Z dt



J ¼ 1;1 J2;1

J1;2 J2;2

 ð4Þ

The eigenvalues of a 2ffi  2 matrix are given by  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l ¼ ð1=2Þ tr + tr2  4det , where tr is the trace (i.e. sum of elements on the main diagonal) and det is the determinant (Strang, 1988). Thus, for Re(l ) , 0, the necessary conditions for stability are (Strang, 1988) tr ¼ J1;1 þ J2;2 , 0

ð5Þ

det ¼ J1;1 J2;2  J2;1 J1;2 . 0

ð6Þ

In later sections, the influence of the functional response on these stability conditions is examined using PZ models based on simplified versions of the FWF

1217

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

NPZ model (“The original NPZ model of FWF” section).

R E S U LT S Analysis of the original FWF model In FWF, two functional responses were contrasted. The first was Ivlev [Ivlev, 1955, equation (7), Fig. 1],   I ¼ 1  eaP Imax

ð7Þ

where the specific ingestion rate increases with phytoplankton density to a maximum value, Imax. The rate of increase is dictated by the parameter a, such that higher a corresponds to a lower prey density at which ingestion is essentially satiated. When parameterized as in (FWF; Table I), ingestion is 50% of Imax when phytoplankton density is 0.7 mmol N m23, and 90% when phytoplankton density is 2.3 mmol N m23 (Fig. 1). These percentages are independent of Imax, and the characteristic shape, a concave downward curve, is also independent of parameter values. The second function considered by FWF was Mayzaud – Poulet [Mayzaud and Poulet, 1978; equation (8), Fig. 1],   I ¼ 1  ebP cP

ð8Þ

which does not satiate, representing ingestion as always increasing with prey density. It was intended to describe

Fig. 1. Functional responses for FWF model. Ingestion rate (day21) as a function of phytoplankton density (mmol N m23) as used in FWF for parameter values given in Table I.

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

observations of copepod ingestion that changed proportionately with phytoplankton at high density, which was attributed to acclimation of their grazing response to increased food availability (Mayzaud and Poulet, 1978). The parameter b affects the response when prey density is low, whereas at high prey density, the response is essentially linear according to the parameter c. When parameterized as in (FWF, Table I), Mayzaud –Poulet ingestion is equal to Ivlev when prey density is 1.0 mmol N m23 (Fig. 1). Below that density, Ivlev yields a higher ingestion rate; above it Mayzaud – Poulet ingestion is higher. Different values for b and c affect the intersection point without altering the characteristic shape of the Mayzaud – Poulet response as a concave upward curve that is nearly linear when phytoplankton density is high. The effects of the two functional responses on the modelled dynamics were contrasted in (FWF, their Fig. 7) and are re-created here (Fig. 2, parameter values in Table I). The most obvious difference, as observed in FWF, is that the Ivlev model results in an effective limit cycle, with predator, prey and nutrients exhibiting large oscillations for months (Fig. 2A), whereas the Mayzaud – Poulet model is highly damped with variables reaching a steady state in about a week (Fig. 2B). Furthermore, the mean values predicted using the Mayzaud – Poulet formulation are notably different than those with Ivlev: phytoplankton density is doubled, zooplankton density increases by almost 50%, nutrient density is halved and both primary production and re-mineralization rates are increased by .50%. A broad range of alternative parameter values result in even more dramatic differences (Franks et al., 1986a and this study). For example, the amplitude of the Ivlev oscillations is increased when (i) phytoplankton growth rate is reduced to account for light effects (e.g. Edwards et al., 2000a, b; Spitz et al., 2003), (ii) total nitrogen is increased to simulate different regions (e.g. Franks and Chen, 1996, 2001) or (iii) mortality rate is decreased (e.g. Edwards et al., 2000a, b; Baird et al., 2006). While there are some parameter sets for which the Ivlev oscillations are damped, the damping effect is fairly sensitive to changes in any single parameter, such that the oscillations tend to return when environmental conditions vary (e.g. Edwards et al., 2000a, this study). Thus, the two functional responses predict fundamentally different dynamics for wide ranges of ecosystems.

Does satiation cause the differing model results? It has long been observed that satiating responses can destabilize ecosystems (May, 1973), and this has been

1218

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Fig. 2. Model results for functional responses used in FWF. Time series of modelled density (mmol N m23) of nutrients (N), phytoplankton (P) and zooplankton (Z) for functional responses shown in Fig. 1. (A) Ivlev and (B) Mayzaud–Poulet.

the stated explanation for the model differences shown in Fig. 2 (Franks et al., 1986a; Davidson, 1996; Franks, 2002). However, as discussed in Analysis of the original the FWF model section, the two functional responses also differ when phytoplankton density is low. The range of phytoplankton densities in the simulations is 0 –0.5 mmol N m23 (Fig. 2), which suggests that the results of oscillating versus constant densities arise from differences in the food-limited aspects of the functional responses. Food-limited functional responses are usually characterized by Type (Holling, 1959). Type 1 ingestion increases linearly with food. I ¼ dP

ð9Þ

such that d is the constant rate of change of ingestion (Fig. 3A). Type 2 responses, which include Ivlev and Michaelis– Menten (also known as Holling Disk, Gentleman et al., 2003), i.e. I¼

P Imax kP þ P

ð10Þ

where Imax is the maximum ingestion rate and kP is the half-saturation constant as in equation (1) for nutrient uptake, exhibit a decreasing rate of change of ingestion as food becomes more abundant such that the relationship is a concave downward curve (Fig. 3B). Type 3 responses are concave downward only once food surpasses a certain level, an inflection point. Below that point, classic Sigmoidal Type 3 responses, which can be expressed in terms of Michaelis– Menten parameters

(Gentleman et al., 2003) as I¼

kP2

P2 Imax þ P2

ð11Þ

exhibit increasing rates of change of ingestion as food becomes more abundant (Fig. 3C), whereas Threshold Type 3 responses, i.e. I¼

Peff Imax kP þ Peff

ð12Þ

where Peff = max(P2Pthresh,0), exhibit no ingestion (Fig. 3C). The Type 2 Ivlev and Michaelis– Menten responses are not mathematically equivalent in that there are no parameter values for which the two functions overlap. They exhibit the same half-saturation value, i.e. phytoplankton density at which ingestion is half its maximum rate when kP = (ln 2)/a, and for such a parameterization the Michaelis –Menten equation predicts slightly higher ingestion rates than Ivlev when P , kP and vice versa (Fig. 3B). Simulations using the Type 2 Michaelis– Menten response result in oscillations just as Type 2 Ivlev did (Fig. 4B), although these are less damped. The Mayzaud – Poulet response does not correspond to any of the classic Types. When phytoplankton density is high, Mayzaud – Poulet is essentially a Type 1 [i.e. equation (8) at high P: IcP], but when phytoplankton density is low, Mayzaud – Poulet is similar to a Sigmoidal Type 3 in that the rate of change of ingestion increases with phytoplankton. The key difference between a classic Sigmoidal Type 3 response and Mayzaud – Poulet is that Type 3 formulations

1219

JOURNAL OF PLANKTON RESEARCH

j VOLUME

30

1220

j NUMBER

11

j PAGES

1215 – 1231

j 2008

Fig. 3. Classic functional responses and associated clearance rates. Classic functional responses describing specific ingestion rate, I, (day21) versus phytoplankton density (mmol N m23): (A) Type 1: Linear, d = 1.5 day21, (B) Type 2: Michaelis– Menten (solid line) and Ivlev (dashed line), kP= 0.7 mmol N m23, a = 1 (mmol N m23)21, Imax = 1.5 day21, (C) Type 3: Sigmoidal (solid line) and Threshold (dashed line), kP = 0.7 mmol N m23, Imax = 1.5 day21, Pthresh = 0.2 mmol N m23. Clearance rates (I/P) for functional responses given in (A–C) versus phytoplankton density: (D) Type 1, (E) Type 2 and (F) Type 3.

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Fig. 4. Model results using classic functional responses. Time series of modelled density (mmol N m23) of nutrients (N), phytoplankton (P) and zooplankton (Z) for classic functional responses as parameterized in Fig. 3. (A) Type 1, (B) Type 2 Michaelis-Menten and (C) Type 3 Sigmoridal.

exhibit satiation. Yet, the model results using the Type 3 Sigmoidal rapidly converge to a stable steady state just as for Mayzaud – Poulet (Fig. 4C). The results for the Threshold Type 3 response are similar to the Sigmoidal Type 3, and steady states are achieved for a wide range of parameter values for both Type 3 responses. In contrast, use of a non-satiating Type 1 ultimately results in a steady state, but the time taken for the initial oscillations to damp is highly dependent on d and can be quite long. For example, when Type 1 is parameterized with a slope similar to the high phytoplankton region of Mayzaud –Poulet, oscillations remain large even after 2 months (Fig. 4A). Thus, it is not satiation that causes the model differences in Fig. 2, nor does a non-satiating response necessarily confer stability to the dynamics.

Control of food-limited dynamics: theoretical derivation and biological interpretation From the preceding section, it is evident that the difference in the model results in Fig. 2 is due to seemingly subtle differences among the food-limited aspects of the functional response. To understand how food-limitation affects the dynamics, we examine the influence of the Type of functional response on the stability conditions for a PZ system [equations (5) and (6)]. For this PZ model, the phytoplankton equation [equation (1)] is simplified by considering the phytoplankton specific growth rate, m, as constant (e.g. nutrient levels are relatively stable), and grazing as the only loss to phytoplankton (e.g. cell death is negligible) such that the phytoplankton equation becomes: dP ¼ mP  IZ dt

ð13Þ

Using equation (2) for the zooplankton equation results in a generalized predator-prey model that is equivalent to Lotka – Volterra (Lotka, 1925; Volterra, 1926) when the functional response is Type 1. The steady state conditions for this PZ model are mP = IZ [from equation (13)] and aI = mZ [from equation (2), when only non-zero values of zooplankton are considered). The corresponding Jacobian [equation (4)] is 

dI m  dP Z I J¼ dI a dP Z aI  mZ   I dI    Z I ¼ P dIdP a dP Z 0



ð14Þ

For this Jacobian matrix, det (J) ¼ aI (dI/dP)Z, which is always positive [i.e. satisfies equation (6)] except when ingestion is satiated and there is no firstorder influence of the functional response on stability (i.e. I = Imax, dI/dP = 0 and Re(l ) = 0). For non-satiated ingestion, the system will be stable when tr (J) ¼ (T/P 2 dI/dP Z,0) [equation (5)], which occurs whenever dI/dP.I/P. While the stability requirement above has been acknowledged in previous studies (e.g. Oaten and Murdoch, 1975a), we demonstrate here that there is biological significance to the mathematical condition. First, we note that the ratio I/P is equal to the clearance rate, C, which is the volume of water filtered per unit time by the zooplankton (Frost, 1972). This allows the stability requirement to be rewritten as dI C .0 dP

ð15Þ

Second, we derive the rate of change of clearance rate with phytoplankton density (i.e. the slope of the

1221

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

clearance rate versus phytoplankton curve):

dC d I ððdI=dPÞ  CÞ ¼ ¼ dP dP P P

ð16Þ

As the term in parentheses in equation (16) is identical to the stability requirement [equation (15)], the sign of dC/dP is identical to the sign of the eigenvalues and therefore has equivalent effect on the dynamical stability. The system is stable when dC/dP . 0 because zooplankton will be able to suppress blooms when clearance rates increase as phytoplankton increases, and because grazing pressure will be reduced when clearance rates decrease as phytoplankton decreases. The system is unstable when dC/dP , 0 because grazing cannot keep phytoplankton growth in check when clearance rates decrease as phytoplankton increases, and because grazing will exacerbate phytoplankton losses when clearance rates increase as phytoplankton decrease. The different Types of food-limited responses exert different influences on the dynamics due to the density-dependence of their associated clearance rates. Type 1 has constant clearance rates (i.e. dC/dP = 0, Fig. 3D), and therefore no first-order influence on stability. Type 2 clearance rates always decrease as phytoplankton increases (Fig. 3E), and are therefore unconditionally unstable. In contrast, Type 3 clearance rates increase with phytoplankton when food is low, and then reach a maximum value after which the clearance rate decreases with further increase in phytoplankton (Fig. 3F). Thus, only Type 3 responses can ever satisfy the stability requirement, and even then, only when the steady-state phytoplankton density is within the low range where clearance rates are increasing. The destabilizing influence of a Type 2 response can be mitigated by other processes. For example, if the zooplankton equation of the PZ model is changed from linear to quadratic closure, equation (2) is replaced by dZ ¼ aIZ  mZ Z 2 dt

ð17Þ

This modification results in corresponding changes to the steady-state values (i.e. aI = mZZ from equation (17), when only non-zero values of zooplankton are considered) and the Jacobian (i.e. J2,2 = aI22mZZ=2mZZ). The associated stability conditions become: dC/dP þ mZ/P.0 [from equation (5)] and dC/dPþ 1/P dI/ dP.0 [from equation (6)]. Thus, quadratic closure provides a mechanism by which the system may be stable

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

in situations for which it is not possible through the functional response alone (i.e. where dC/dP.0). Again, this mathematical finding can be interpreted biologically when the closure term represents predation. Predation rates can be characterized as the product of the predators’ functional response (Ipred) and their biomass (B). Thus, if we consider predators whose biomass is essentially constant over the period of interest, and defining d = mZ/B, then (i) linear closure corresponds to predators with a Type 1 functional response on zooplankton (i.e. Ipred = dZ), whereas (ii) quadratic closure represents a non-standard functional response (i.e. Ipred = dZ 2), which is similar to Mayzaud –Poulet [equation (8)] in that it is like a Sigmoidal Type 3 functional response at low zooplankton density and does not satiate. With this interpretation, the damping effect of quadratic closure is clear: it allows the predators’ clearance rates to increase when zooplankton biomass increases (and vice versa), thereby enabling the predators to rapidly respond to any changes in zooplankton in an analogous manner to what was discussed above for zooplankton feeding on phytoplankton.

The influence of clearance rates in NPZ models: a sensitivity study As shown in the previous section, it is possible for the PZ system to exhibit stability even when dC/dP , 0 depending on the manner in which other processes affect the sign of the Jacobian’s eigenvalues. This is also true for more complex systems, such as the NPZ model of FWF (“The original NPZ model of Franks et al. (1986a) Section”), for which the eigenvalues are affected by nutrient-dependent growth rates and phytoplankton cell death in addition to zooplankton clearance rates. While it is possible to conduct stability analyses on more complex models (e.g. Busenberg et al., 1990; Newberger et al., 2003), determination of the conditions for negative eigenvalues is not as straightforward as above (i.e. Jacobian terms must satisfy Routh – Hurwitz criteria for larger systems: Nisbet and Gurney, 1982) and analyses are still limited to the asymptotic behaviour of first-order effects. In order to examine the full effect of nonlinearities and demonstrate the importance of clearance rates we use sensitivity studies where parameters and functions of the original FWF NPZ model (Table I) are varied and the time-series of results are examined. The simulations in the previous sections have already shown that the density-dependence of clearance rates strongly influences the dynamics in more complex systems. For example, it was shown that the ecosystem is unstable with Type 2 functional responses where dC/ dP , 0 and stable with Type 3 where dC/dP . 0

1222

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

(Fig. 4B versus C). Clearance rates also explain why the Michaelis– Menten Type 2 response generated larger amplitude oscillations than the Ivlev Type 2 (Fig. 4B versus Fig. 2B), as the slope of clearance rate versus phytoplankton curve is much flatter for Ivlev (i.e. dC/ dP is more negative for Michaelis– Menten, Fig. 3E). Furthermore, they explain the stabilizing effect of Mayzaud – Poulet at low phytoplankton (i.e. dC/dP = cb e 2bP . 0; Fig. 2A) and why use of the non-satiating Type 1 resulted in a weakly stable system (i.e. ecosystem stability is due to other processes as Type 1 has no firstorder influence, Fig. 4A). The simulations discussed above were all generated using the original parameter set (Table I). There are two ways that use of different parameter values can affect the system’s stability. First, changing parameters alter the steady-state value of phytoplankton, which alters the point on the functional response at which dC/dP is evaluated. Second, if the parameter being changed relates to the functional response, it will additionally alter the shape of the functional response, such that dC/dP is altered for any particular phytoplankton density. As mentioned earlier, varying certain ecological parameters (e.g. decreasing mmax or mZ, or increasing total mass) amplifies the oscillations resulting from use of a Type 2 functional response. In this model, all these scenarios result in a lower steady-state value of phytoplankton, which results in a more negative value of dC/ dP (Fig. 5A and B), and therefore a stronger destabilizing influence of the Type 2 functional response. Variation in functional response parameter values can also amplify the inherent instability of the Type 2 response. For example, when either kP is halved or Imax doubled, the oscillations grow so large that plankton go extinct. The reason for such a dramatic effect is that these parameter changes lead to both a lower steady-state phytoplankton density and a higher nonlinearity of the functional response such that dC/dP at steady state is significantly more negative (Fig. 5A and B). Extinction is not seen with a Type 3 response: even when kP decreases or Imax increases by an order of magnitude, oscillations never develop and plankton persist. The reason for the difference from the Type 2 simulation is that such variation in parameters with Type 3 responses make dC/dP more strongly positive (i.e. more stabilizing) at lower phytoplankton density (Fig. 5C and D). It is possible to generate stable ecosystem dynamics with a Type 2 functional response by changing parameters. For example, if either kP is doubled or Imax is halved, oscillations are damped and the model results in constant densities after several months. Essentially, both

of these scenarios flatten the functional responses and increase the steady-state phytoplankton values such that dC/dP at steady state is less negative than in the original case (Fig. 5A and B). Therefore, the destabilizing influence of the Type 2 functional response is reduced sufficiently that it can be overcome by the stabilizing effect of other modelled processes. As with the PZ model (see above), use of quadratic closure may stabilize models even when dC/dP , 0, but the effect is dependent on specific parameter values. For example, using quadratic closure (mortality = mZZ 2) instead of linear closure in the FWF, NPZ model (Table I) stabilizes the ecosystem even when a Type 2 response is used, but only after several months of simulation. The parameter mZ must be increased by more than a factor of 3 for the damping to be comparable to that seen when Type 3 grazing and linear closure are employed (Fig. 4C). Thus, while closure is important, it does not always remove oscillations entirely and cannot completely negate the influence of functional response characterization.

Effect of acclimation time lags Type may not be the only aspect of the functional response that can affect ecological stability, as time lags can be destabilizing (May, 1973; Steele and Henderson, 1981; Ruan, 2001) and copepods require an acclimation period on the order of 1 –6 days to adapt their feeding behaviour to changes in food conditions (Mayzaud and Poulet, 1978; Davis and Alatalo, 1992; Hirche et al., 1997, Niehoff, 2004). FWF modified their functional response to explore how such time lags could affect modelled dynamics and noted only a small destabilizing effect. However, their approach allowed for a large instantaneous adjustment in ingestion rates that may be unrealistic (see Appendix). Here, an acclimation time lag was implemented in the FWF NPZ model (Table I) by replacing P(t) in the functional response with Pavg,t (t), a running average of phytoplankton density over the past t days. Thus, when phytoplankton density changes, zooplankton modify their ingestion rates gradually, only reaching the corresponding functional response value after t days (Fig. A1). Because Pavg,t (t) may be higher than P(t), the model also included a conditional statement that limited predicted consumption to be no greater than the amount of phytoplankton available. To stress the possible destabilizing effect of acclimation lags, the model was set up for the most stable conditions: i.e. Type 3 functional response and quadratic closure. Initial variable values were set to their steady states, and Pavg,t (0) was set to the steady-state value for phytoplankton. For t , 1.5 days, the system maintains

1223

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

Fig. 5. Effect of parameters on density-dependence of clearance rates. Variation of rate of change of clearance rate with phytoplankton (dC/dP, i.e. slope of C versus P curve) for Michaelis–Menten Type 2 (Fig. 3E) and Sigmoidal Type 3 (Fig. 3F) with parameters (A and C) kP, and (B and D) Imax at original, half and double values.

its steady state. However, for longer lags, oscillations develop and grow such that for t . 2 days zooplankton drive the phytoplankton to extinction within 2 – 3 weeks (Fig. 6). The destabilizing influence of acclimation lags is robust to changes in the functional response parameters that increase dC/dP at steady state (e.g. kP is halved or Imax is doubled, Fig. 5C and D). In contrast, when parameter values are changed such that dC/dP is reduced at steady state, acclimation lags become less important. For example, when kP is doubled or Imax is halved (Fig. 5C and D), oscillations do not develop until t . 5 days, and these are sufficiently damped that extinction does not occur within 2 months even for t as large as 10 days. The destabilizing effect and parameter sensitivity of acclimation lags is made clear by considering the corollary effect on instantaneous clearance rates. Using chain rule, the density dependence of clearance rates can be

re-written for responses that are in terms of Pavg,t, i.e. dC dC dPavg;t ¼ dP dPavg;t dP

ð18Þ

It can be seen from equation (18) that even when phytoplankton density is in the range where the Type 3 functional response is stable (dC/dPavg,t . 0), clearance rates will exert a destabilizing influence (i.e. dC/dP , 0) whenever dPavg,t/dP , 0. This latter condition will occur whenever the present and historical phytoplankton temporal trends are out of phase (i.e. opposite signs of dP/dt and dPavg,t/dt), for example when phytoplankton density increases more slowly or decreases more rapidly than it did t days prior. Thus, the critical lag observed in our simulations is due to the fact that t must reach the inherent time scale of prey variability

1224

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Fig. 6. Effect of acclimation time lag. Time series of modelled density (mmol N m23) of nutrients (N), phytoplankton (P) and zooplankton (Z) for FWF model with Type 3: Sigmoidal functional response with Pavg,t (t) where t = 3 days, quadratic closure and initial conditions based on steady-state conditions: N0 = 0.76 mmol N m23, P0 = 0.32 mmol N m23 and Z0 = 0.92 mmol N m23.

before the density dependence of the clearance rates becomes negative. The parameter sensitivity results from the fact that jdC/dPj increases when dC/dPavg,t increases (Fig. 5C and D) and therefore the functional response exerts a stronger destabilizing influence when zooplankton are more severely food-limited. However, when dC/dPavg,t is sufficiently small, other processes may be able to stabilize the system even when dPavg,t/ dP , 0.

DISCUSSION Here we have shown that the functional response plays a pivotal role in models of planktonic ecosystems. We demonstrated that variations in its mathematical form significantly affect results both qualitatively (oscillating versus constant densities) and quantitatively (mean biomass and production), including extinction versus persistence. In fact, changes in the food-limited and/or transient nature of the functional response arguably have a more profound effect on ecological stability than changes in mortality terms (closure). Thus, the choice of response can fundamentally affect ecological interpretations and model predictions. Our analyses revealed that the influence of the functional response can be explained by the density dependence of the clearance rates, C: grazing acts to suppress blooms and affords prey refuges when dC/dP . 0, whereas it acts to exacerbate changes in phytoplankton

when dC/dP , 0. Clearance rates, which are related to feeding preferences (Gentleman et al., 2003), are affected by the zooplankton’s perceptual abilities (chemical, visual and mechanical) and feeding modes (i.e. suspension versus ambush), as well as turbulence (Saiz and Kiorboe, 1995; Kiorboe et al., 1996). In order to make appropriate choices, modellers need to understand how different feeding behaviours alter zooplankton clearance rates. Common practices suggest the role of clearance rates and resulting dynamical effects of satiation, food-limitation and acclimation are not appreciated, which may limit the utility of current ecosystem models. While our finding is technically consistent with historical views that satiation is destabilizing (i.e. when I = Imax, C = Imax/P and therefore dC/dP , 0), we showed this effect is relatively weak even at moderate phytoplankton densities (i.e. jdC/dPj is small for P . 2kP , Figs. 3 and 5). Furthermore, we refuted historical views about the stabilizing effect of non-satiating responses by (i) demonstrating that it is not the nonsatiating nature of the Mayzaud – Poulet response that is responsible for its stabilizing influence in FWF, and (ii) showing that at high prey density Mayzaud –Poulet is like a Type 1, which has no first-order effect on stability (“Control of food-limited dynamics: theoretical derivation and biological interpretation” section). Additionally, use of non-satiating responses may result in unrealistic ingestion rates at high prey density. For example, at the phytoplankton bloom densities encountered by the copepods on which Mayzaud – Poulet based their formulation, and using parameters in Table I, modelled ingestion is 15 day21—roughly an order of magnitude higher than observed maximum rates (Frost, 1972; Hansen et al., 1997). Similarly, predation rates in many models may be greatly overestimated when zooplankton density is high, as both linear and quadratic closure can be construed as non-satiating responses for the predators (“Control of food-limited dynamics: theoretical derivation and biological interpretation” section). We have shown that the clearance rates associated with different food-limited Types of functional responses have fundamentally different density-dependence (Fig. 3) and thereby exert fundamentally different influences on ecological dynamics. This explains the findings of previous studies focusing on functional response and stability. For example, the “interspecific interaction” term in the classic Lotka – Volterra predator-prey model is equivalent to a constant clearance rate (i.e. ingestion is Type 1, dC/dP = 0) and this model is well-known to exhibit neutral stability at the non-trivial steady state (e.g. May, 1973). Furthermore, the destabilizing

1225

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

influence of the Type 2 response has been shown to be sensitive to both parameters and functional form (e.g. Myerscough et al., 1996; Fussman and Blasius, 2005), which can be understood in terms of the how these parameters and/or functions affect the extent to which dC/dP is negative (“The influence of clearance rates in NPZ models: a sensitivity study” section). The ability of Type 3 functions to stabilize is also known to be dependent on parameters (Steele, 1974; Myerscough et al., 1996), which can be understood in terms of how the parameters affect whether and to what extent dC/dP is positive at steady state (Fig. 5C and D). The differing dynamic influence of the Types also explains why quadratic closure is generally more stabilizing than linear closure in that at low zooplankton density, these closure schemes are akin to Type 1 (dCpred/dZ = 0) versus Type 3 (dCpred/dZ . 0) responses for the predators (“Control of food-limited dynamics: theoretical derivation and biological interpretation” section). Despite the fact that the destabilizing influence of the Type 2 functional response is sensitive to parameters, modellers seldom provide biological justification for their choice of a for the Ivlev functional response. This may be because empiricists tend to use Michaelis – Menten equations to fit their data (e.g. Hansen et al., 1997; Hirst and Bunker, 2003), and modellers may not recognize the relationship between a and kP. Commonly used Ivlev constants of 0.2– 0.3 (mmol N m23)21(Franks and Chen, 1996, 2001; Franks and Walstad, 1997; Edwards et al., 2000a, b; Baird et al., 2006) correspond to kP 200 mg C m23 (using Redfield C:N), which is at the upper end of the range in measured values for copepods (e.g. 30– 300 mg C m23: Hansen et al., 1997; Hirst and Bunker, 2003 using C:Chl=50). Some applications have even used a = 0.06 (mmol N m23)21 (Newberger et al., 2003; Spitz et al., 2003), for which ingestion is modelled as severely food-limited at bloom-level densities (i.e..500 mg C m23). Use of an unrealistically high kP (or low a) flattens the Type 2 functional response thereby reducing the destabilizing nature of the clearance rates (Fig. 5) and underestimating the influence of grazing on the system dynamics. Even when Type 2 functional response parameter values are realistic, they can vary with both physical and biological factors. For example, maximum grazing rates (Imax) for microzooplankton are typically 2 – 10 times higher than mesozooplankton (Hansen et al., 1997). Imax also varies with temperature (i.e. Q10 = 2 – 4: Campbell et al., 2001; Hirst and Bunker, 2003), and can depend on food quality and size (Frost, 1972). Models often describe kP (or a) as increasing with the density of alternative prey (e.g. Multiple-prey Michaelis– Menten: Gentleman et al., 2003), predators (DeAngelis et al.,

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

1975) and turbulence (Rothschild and Osborne, 1988). Because jdC/dPj is greater for higher maximum grazing rates (Imax) and lower half-saturation constants (kP; Fig. 5), modelled Type 2 responses accounting for these factors will result in ecosystem dynamics being more prone to oscillations for smaller taxa, warmer temperatures, more nutritious food, more selective diets, decreased competition and/or stratified waters (e.g. Edwards et al., 2000a, this study). Note that the same environmental variations are unlikely to induce oscillations in models using a Type 3 functional response for zooplankton, as the stabilizing influence of such responses is robust to order of magnitude changes in Imax and kP. Thus, the fundamental difference in food-limited Type 2 versus Type 3 responses may lead to contrasting conclusions about the importance of physical and biological factors affecting the functional response parameters. Furthermore, when models are “tuned” for stability, models using Type 2 will necessarily result in a higher relative importance of other stabilizing processes such as predation (Steele and Henderson, 1981, 1992, 1995) or mixing (Edwards et al., 2000b). Results may also be biased by the way models represent zooplankton acclimation to changes in food. For example, the Mayzaud –Poulet formulation is often referred to as an acclimation response (e.g. Wroblewski and Richman, 1987; Doney et al., 1996; McClain et al., 1996; Leonard et al., 1999; Botte and Kay, 2000; Friedrichs and Hofmann, 2001; Franks, 2002; Tian 2006), thereby implying that acclimation provides a stabilizing influence on ecosystem dynamics. However, this response was intended to characterize behavioural changes over time scales longer than those associated with prey variability (Mayzaud and Poulet, 1978). Thus, the stabilizing aspect of the Mayzaud –Poulet response (i.e. the fact that it characterizes feeding for which dC/dP . 0) does not account for transient behaviour in dynamic prey fields. Furthermore, the unrealistic implementation of such behaviours in FWF led to the false conclusion that time lag effects are negligible. We demonstrated that acclimation lags of 2 days could destabilize a model with Type 3 functional response and quadratic closure due to the effect on the zooplankton clearance rates (“Effect of acclimation lags” section). We showed that the influence of acclimation lags depends on the magnitude of two factors: the relative scale of historical versus present prey variation (i.e. dPavg,t/dP) and the density-dependence of the clearance rates (“Effect of acclimation lags” section). The latter is related to the extent to which zooplankton are food-limited, and explains why we found the destabilizing effect of lags was robust to both decreases in the

1226

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

half-saturation constant (kP) and increases in maximum ingestion rates (Imax). This suggests acclimation lags could be dynamically important for food-limited mesozooplankton in a broad range of environments, although the fact that the oscillations were large enough to drive plankton to extinction indicates processes not modelled must contribute to the maintenance of their prey populations. While we also found parameter values for which acclimation lags were weak and could be overcome by other processes, such ecological stability would not be sustainable if other factors caused phytoplankton variability to increase and thereby magnify the destabilizing influence of clearance rates [i.e. dC/dP will be more negative when dPavg,t/dP larger, equation (17)]. Acclimation lags may not be as significant for microzooplankton as they are for mesozooplankton. Although microzooplankton have higher Imax (i.e. tendency to destabilize), they are typically modelled with higher kP (tendency to stabilize) and likely exhibit shorter time periods before behaviour adjusts to the new conditions. Ultimately, the importance of including acclimation lags for zooplankton in any particular environment depends on the feeding behaviour (i.e. Type) and whether the time scale for behavioural adjustments is larger than that of prey variability. Thus, the choice of Type and transient nature of the functional response is critical, as it can affect conclusions about the role of grazing, bottom-up versus top-down control and the influence of climate. However, the appropriate choice of functional response is not always apparent. There are often insufficient data on which to base decisions due to the challenges of conducting experiments at low food levels, such as maintaining uniform prey density in the face of consumption, increased mortality due to poor condition of the animals, etc. Uncertainty in the data often enables reasonable fits with different Types and there may be no statistical basis for choosing one Type over another (Mullin, 1975). Many empiricists tend to use a standard Type 2 function without exploring other forms, although there are numerous biological reasons for Type 3 responses including switching between alternative prey types (e.g. Hassell et al., 1977; Strom et al., 2000; Gentleman et al., 2003). Furthermore, since zooplankton are typically acclimated to food levels for a minimum of 24 h in functional response experiments (e.g. Frost, 1972; Verity, 1991), these empirical relationships provide no information on the transient response during the period of acclimation. The practice of using formulations which are more stabilizing, as is often the stated justification for choosing Mayzaud – Poulet

or Type 3 responses, may also be questionable. Oscillations do occur in the field (e.g. McCauley and Murdoch, 1987), and even when observations indicate stable dynamics, that stability may derive from other factors, such as density-dependent zooplankton mortality (this study, “The influence of clearance rates in NPZ models: a sensitivity study” section), mixing (Edwards et al., 2000b) and/or feeding on additional prey types (e.g. Armstrong, 1999; Lima et al., 2002). What should modellers do to ensure the functional response formulation does not limit the utility of their ecosystem model? As they may overestimate ingestion, we recommend against use of unrealistic non-satiating responses for zooplankton (e.g. Mayzaud – Poulet or Type 1) or for their predators (e.g. linear or quadratic closure). Instead, we advocate use of satiating responses for zooplankton (e.g. Type 2 or Type 3) and their predators (Steele and Henderson, 1992, Neubert et al., 2004). As parameter values affect clearance rates and therefore ecological stability, functional response parameters should not simply be recycled from previous models. Instead, values should be based on data that is characteristic of the taxa and environment being studied, which may necessitate separate responses for microzooplankton and mesozooplankton and/or multiple prey types. In addition, sensitivity studies should not be restricted to traditional analyses that only vary parameters over limited ranges (e.g. Fasham, 1995; Evans, 1999, Wainwright et al., 2007), as the dynamical influence of particular parameters depends on Type, and variation in model structure (e.g. Type, time delays, etc.) can affect results at least as much as parameter values. To more accurately quantify uncertainty, modellers should conduct sensitivity studies that vary both parameters and functional forms (e.g. Leising et al., 2003). Our simple model of acclimation is not mechanistically based, and we suggest modellers develop formulations that explain behavioural changes in more realistic ways, including consideration of internal nutrient reserves separate from ambient prey levels akin to what is done with quota models of phytoplankton (Jeschke et al., 2002; Flynn, 2003; Mitra et al., 2007; Mitra and Flynn, 2007).

FUNDING This work was supported by grants from the U.S. National Science Foundation (Biological Oceanography—0222309) and Natural Sciences and Engineering Research Council of Canada and is GLOBEC contribution 599.

1227

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

AC K N OW L E D G E M E N T S We would like to thank Kevin Flynn, Charles Hannah, Moritz Lehmann, Chrissy Galloway, Romero Advincula and three anonymous reviewers for their encouraging discussions and editorial feedback on earlier versions of this manuscript.

APPENDIX FWF implemented a time delay in the Mayzaud – Poulet functional response [equation (8)] by replacing the “cP” term with a dynamic variable A: I ¼ Að1  ebP Þ

ðA1Þ

and solved for A based on its rate of change dictated by time lag t dA cP  A ¼ t dt

ðA2Þ

where A(0) = cP(0). When t less than or equal to the numerical integration time step, this approach is identical to the use of the Mayzaud – Poulet ingestion rate without a time delay. For larger t, there is a time lag before A adjusts to its value at the new phytoplankton density. As the time lag only affects A, the exponential term in equation (A1) still changes instantaneously with

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

phytoplankton. This results in a modelled time delay of ingestion that is not intuitive. To illustrate this, consider the simple scenario where phytoplankton density has been constant at Pold for a length of time . t, and then is suddenly changed to a new value, Pnew, where it remains, as might occur in a laboratory experiment. The acclimation transient response therefore models the change from I(Pold) to I(Pnew). Use of the FWF formulae results in modelled ingestion rates that undergo a large instantaneous change equalling 30% of DI. Subsequently, the modelled ingestion approaches I(Pnew) at a decaying exponential rate. The decay is slow such that the modelled change in ingestion is only 75% of DI at t and only 91% at 2t. In contrast, using our time-delay approach (i.e. replacing P(t) with Pavg,t (t) as in section ‘Effect of acclimation time lags’) for the Mayzaud – Poulet functional response results in a linear interpolation between the two predicted ingestion rates, such that at t the modelled ingestion rate has fully acclimated to the new phytoplankton level (Fig. A1). Note that, when t = 3 days, this approach to acclimation results in ingestion taking 1 day to attain the same adjustment that the FWF model does instantaneously, and the predicted rate of ingestion is less than FWF for the first 2 days.

REFERENCES Arhonditsis, G. B., Adams-Vanharn, B. A., Neilsen, L. et al. (2006) Evaluation of the current state of mechanistic aquatic biogeochemical modeling: citation analysis and future perspectives. Environ. Sci. Technol., 40, 6547–6554. Armstrong, R. A. (1999) Stable model structures for representing biogeochemical diversity and size spectra in plankton communities. J. Plankton Res., 21, 445–464. Aydin, K. Y., McGarlane, G. A., King, J. R. et al. (2005) Linking ocean food webs to coastal production and growth rates of Pacific salmon (Oncorhynchus spp.), using models on three scales. Deep-Sea Res. II, 52, 757–780. Batchelder, H. P., Edwards, C. A. and Powell, T. M. (2003) Individual-based models of copepod populations in coastal upwelling regions: implications of physiologically and environmentally influenced diel vertical migration on demographic success and nearshore retention. Prog. Ocean., 53, 307– 333. Baird, M. E., Timko, P. G., Suthers, I. M. et al. (2006) Coupled physical– biological modelling study of the East Australian Current with idealised wind forcing. Part I: biological model intercomparison. J. Mar. Syst., 59, 249 –270.

Fig. A1. Models of acclimation time delays. Modelled ingestion rate (day21) versus time for case where P = Pold = 0.5 mmol N m23 until t = 0, when it changes to Pnew = 1.0 mmol N m23, and ingestion rates are determined by Mayzaud–Poulet functional responses with instantaneous response to change in prey (i.e. no acclimation) and acclimation modelled as in FWF and in this study.

Botte, V. and Kay, A. (2000) A numerical study of plankton population dynamics in a deep lake during the passage of the Spring thermal bar. J. Mar. Syst., 26, 367– 386. Busenberg, S., Kumar, S. K., Austin, P. et al. (1990) The dynamics of a model of a plankton-nutrient interaction. Bull. Math. Biol., 52, 677 –696.

1228

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Campbell, R. G., Wagner, M. M., Teegarden, G. J. et al. (2001) Growth and development rates of the copepod Calanus finmarchicus reared in the laboratory. Mar. Ecol. Prog. Ser., 221, 161– 183.

Franks, P. J. S., Wroblewsi, J. S. and Flierl, G. R. (1986b) Prediction of phytoplankton growth in response to the frictional decay of a warmcore ring. J. Geophys. Res. Oceans 91, 7603– 7610.

Carlotti, F. and Wolf, U. (1998) Population dynamics model of copepod species coupled with a 1D model of phytoplankton dynamics. Fish. Oceanogr., 7, 191– 204.

Franks, P. J. S. and Chen, C. (1996) Plankton production in tidal fronts; a model of Georges Bank in summer. J. Mar. Res., 54, 631–651.

Christian, J. R., Verschell, M. A., Murtugudde, R. et al. (2002) Biogeochemical modelling of the tropical Pacific Ocean. I: seasonal and interannual variability. Deep-sea res. II Top. Stud. Oceanogr., 49, 509–543. Davidson, K. (1996) Modelling microbial food webs. Mar. Ecol. Prog. Ser., 145, 279 –296. Davis, C. S. and Alatalo, P. (1992) Effects of constant and intermittent food supply on life-history parameters in a marine copepod. Limnol. Oceanogr., 37, 1618– 1639. DeAngelis, D. L., Goldstein, R. A. and O’Neill, R. V. (1975) A model for trophic interaction. Ecology, 56, 881 –892. Denman, K. L. and Pen˜a, M. A. (1999) A coupled 1-D biological/ physical model of the northeast subarctic Pacific Ocean with iron limitation. Deep-Sea Res. II, 46, 2877–2908. Denman, K. L. and Pen˜a, M. A. (2002) The response of two coupled 1-D mixed layer/planktonic ecosystem models to climate change in the NE Subarctic Pacific Ocean. Deep-Sea Res. II, 49, 5739– 5757. Doney, S. C., Glover, D. M. and Najjar, R. G. (1996) A new coupled, one-dimensional biological-physical model for the upper ocean: Applications to the JGOFS Bermuda Atlantic Time-series Study (BATS) site. Deep Sea Res. II, 43, 591– 624.

Franks, P. J. S. and Chen, C. (2001) A 3-D prognostic numerical model study of the Georges Bank ecosystem. Part II: biologicalphysical model. Deep-Sea Res. II, 48, 457– 482. Franks, P. J. S. and Walstad, L. J. (1997) Phytoplankton patches at fronts: a model of formation and response to wind events. J. Mar. Res., 55, 1–29. Friedrichs, M. A. M. and Hofmann, E. E. (2001) Physical control of biological processes in the central equatorial Pacific Ocean. Deep-Sea Res. I, 48, 1023–1069. Frost, B. W. (1972) Effects of size and concentration of food particles on the feeding behaviour of the marine planktonic copepod Calanus pacificus. Limnol. Oceanogr., 17, 805–815. Frost, B. W. (1987) Grazing control of phytoplankton stock in the open subarctic Pacific Ocean: a model assessing the role of mesozooplankton, particularly the large calanoid copepods Neocalanus spp. Mar. Ecol. Prog. Ser., 39, 49– 68. Frost, B. W. and Franzen, N. C. (1992) Grazing and iron limitation in the control of phytoplankton stock and nutrient concentration: a chemostat analogue of the Pacific equatorial upwelling zone. Mar. Ecol. Progr. Ser., 83, 291– 303. Fussman, G. F. and Blasius, B. (2005) Community response to enrichment is highly sensitive to model structure. Biol. Lett., 1, 9– 12.

Edwards, A. M. and Yool, A. (2000) The role of higher predation in planktonic population models. J. Plankton Res., 22, 1085– 1112.

Gentleman, W., Leising, A., Frost, B. et al. (2003) Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics. Deep-Sea Res. II, 50, 2847– 2875.

Edwards, C. A., Batchelder, H. P. and Powell, T. M. (2000a) Modeling microzooplankton and macrozooplankton dynamics within a coastal upwelling system. J. Plankton Res., 22, 1619–1648.

Hansen, J., Bjornsen, P. K. and Hansen, B. W. (1997) Zooplankton grazing and growth: scaling within the 2 –2000 mm body size range. Limnol. Oceanogr., 42, 687– 704.

Edwards, C. A., Powell, T. M. and Batchelder, H. P. (2000b) The Stability of a NPZ model subject to realistic levels of ocean mixing. J. Mar. Res., 58, 37–60.

Hassell, M. P., Lawton, J. H. and Beddington, J. R. (1977) Sigmoid functional responses and invertebrate predators and parasitoids. J. Anim. Ecol., 46, 249– 262.

Evans, G. T. (1977) Functional response and stability. Am. Nat., 111, 799–802.

Hinckley, S., Napp, J. M. and Hermann, A. J. (2007) Simulation of physically-mediated variability in prey resources for a larval fish—a 3-dimensional NPZ model. Fish. Oceanogr. in press.

Evans, G. T. (1999) The role of local models and data sets in the Joint Global Ocean Flux Study. Deep-Sea Res. I, 46, 1369–1389. Evans, G. T. and Parslow, J. S. (1985) A model of annual plankton cycles. Biol. Oceanogr., 3, 327– 347.

Hirche, H.-J., Meyer, U. and Niehoff, B. (1997) Egg production of Calanus finmarchicus: effect of temperature, food and season. Mar. Biol., 127, 609 –620.

Fasham, M. J. R. (1995) Variations in the seasonal cycle of biological production in subarctic oceans: a model sensitivity analysis. Deep-Sea Res. I, 42, 1111–1149.

Hirst, A. G. and Bunker, A. J. (2003) Growth of marine planktonic copepods: global rates and patters in relation to Chlorophyll a, temperature and body weight. Limnol. Oceanogr., 48, 1988– 2010.

Fasham, M. J. R., Ducklow, H. W. and McKelvie, S. M. (1990) A nitrogen-based model of plankton dynamics in the oceanic mixed layer. J. Mar. Res., 48, 591– 639.

Holling, C. S. (1959) The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Canad. Entomol., 91, 293– 320.

Flynn, K. (2003) Modelling multi-nutrient interactions in phytoplankton; balancing simplicity and realism. Prog. Oceanogr., 56, 249–279.

Holling, C. S. (1965) The functional response of predators to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Can., 45, 5 –60.

Franks, P. J. S. (2002) NPZ models of plankton dynamics: their construction, coupling to physics, and application. J. Oceanogr., 58, 379–387.

Hood, R. R., Kohler, K. E., McCreary, J. P. et al. (2003) A fourdimensional validation of a coupled physical-biological model of the Arabian Sea. Deep-Sea Res. II, 50, 2917–2945.

Franks, P. J. S., Wroblewski, J. S. and Flierl, G. R. (1986a) Behaviour of a simple plankton model with food-level acclimation by herbivores. Mar. Biol., 91, 121–129.

Ivlev, V. S. (1955) Experimental Ecology of the Feeding of Fishes. Pischepromizdat, Moscow, 302 pp. (translated from Russian by D. Scott. New Haven, Yale University Press. 1961).

1229

JOURNAL OF PLANKTON RESEARCH

j

VOLUME

30

j

NUMBER

11

j

PAGES

1215 – 1231

j

2008

Jeschke, J. M., Kopp, M. and Tollrian, R. (2002) Predator functional responses: discriminating between handling and digesting prey. Ecol. Monogr., 72, 95–112.

Neubert, M. G., Klanjscek, T. and Caswell, H. (2004) Reactivity and transient dynamics of predator-prey and food web models. Ecol. Model., 179, 29–38.

Keen, T. R., Kindle, J. C. and Young, D. K. (1997) The interaction of southwest monsoon upwelling, advection and primary production in the northwest Arabian Sea. J. Mar. Syst., 13, 61–82.

Newberger, P. A., Allen, J. S. and Spitz, Y. H. (2003) Analysis and comparison of three ecosystem models. J. Geophys. Res. Oceans, 108, Art no. 3061.

Kiørboe, T., Saiz, E. and Viitasalo, M. (1996) Prey switching behaviour in the planktonic copepod Acartia tonsa. Mar. Ecol. Prog. Ser., 143, 65– 75.

Niehoff, B. (2004) The effect of food limitation on gonad development and egg production of the planktonic copepod Calanus finmarchicus. J. Exp. Mar. Biol. Ecol., 307, 237–259.

Le Fouest, V., Zakardjian, B., Saucier, F. J. et al. (2005) Seasonal versus synoptic variability in planktonic production in a high-latitude marginal sea: The Gulf of St. Lawrence (Canada). J. Geophys. Res. Oceans, 110, Art No. C09012.

Nisbet, R. M. and Gurney, W. S. C. (1982) Modelling Fluctuating Populations. John Wiley and Sons, New York.

Leising, A. W., Gentleman, W. C. and Frost, B. W. (2003) The threshold feeding response of microzooplankton within Pacific highnitrate low-chlorophyll ecosystem models under steady and variable iron input. Deep-Sea Res. II, 50, 2877–2894.

Oaten, A. and Murdoch, W. M. (1975b) Switching, functional response and stability in predator-prey systems. Am. Nat., 109, 299–318.

Leonard, C. L., McClain, C. R., Murtugudde, R. et al. (1999) An iron-based ecosystem model of the central equatorial Pacific. J. Geophys. Res. Oceans, 104, 1325– 1341.

Oaten, A. and Murdoch, W. M. (1975a) Functional response and stability in predator-prey systems. Am. Nat., 109, 289–298.

Rothschild, B. J. and Osborn, T. R. (1988) Small-scale turbulence and plankton contact rates. J. Plankton Res., 10, 465–474. Ruan, S. G. (2001) Oscillations in plankton models with nutrient recycling. J. Theoretical Biol., 208, 15– 26.

Lima, I. D., Olson, D. B. and Doney, S. C. (2002) Intrinsic dynamics and stability properties of size-structured pelagic ecosystem models. J. Plankton Res., 24, 533– 556.

Saiz, E. and Kiørboe, T. (1995) Predatory and suspension feeding of the copepod Acartia tonsa in turbulent environments. Mar. Ecol. Prog. Ser., 122, 147– 158.

Lotka, A. J. (1925) Elements of Physical Biology. Williams & Wilkins Co, Baltimore.

Salihoglu, B. and Hofmann, E. E. (2007) Simulations of phytoplankton species and carbon production in the equatorial Pacific Ocean 1. Model configuration and ecosystem dynamics. J. Mar. Res., 65, 219–273.

May, R. M. (1973) Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton, New Jersey. Mayzaud, P. and Poulet, S. A. (1978) The importance of the time factor in the response of zooplankton to varying concentrations of naturally occurring particulate matter. Limnol. Oceanogr., 23, 1144– 1154. McCauley, E. and Murdoch, W. W. (1987) Cyclic and stable populations: plankton as paradigm. Am. Nat., 129, 97– 121. McClain, C. R., Arrigo, K., Tai, K. S. et al. (1996) Observations and simulations of physical and biological processes at ocean weather station P 1951–1980. J. Geophys. Res. Oceans, 101, 3697–3713. Michaelis, L. and Menten, M. (1913) Invertinwirkung. Biochem. Z., 49, 333–369.

Die

Kinetik

der

Mitra, A. and Flynn, K. J. (2007) Importance of interactions between food quality, quantity and gut transit time on consumer feeding, growth and trophic dynamics. Am. Nat., 169, 632– 646. Mitra, A., Flynn, K. J. and Fasham, M. J. R. (2007) Accounting for grazing dynamics in nitrogen– phytoplankton–zooplankton (NPZ) models. Limnol. Oceanogr., 52, 649–661. Monod, J. (1942) Recherches sur la croissance des cultures bacteriennes. Hermann et Cie, Paris. Mullin, M. M. (1975) Ingestion by planktonic grazers as a function of concentration of food. Limnol. Oceanogr., 20, 259– 262.

Spitz, Y. H., Newberger, P. A. and Allen, J. S. (2003) Ecosystem response to upwelling off the Oregon coast: behaviour of three nitrogen-based models. J. Geophys. Res. Oceans, 108, Art. No. 3062. Steele, J. H. (1974) The Structure of Marine Ecosystems. Harvard University Press, Cambridge, Massachusetts. Steele, J. H. and Henderson, E. W. (1981) A simple plankton model. Am. Nat., 117, 676– 691. Steele, J. H. and Henderson, E. W. (1992) The role of predation in plankton models. J. Plankon Res., 14, 157–172. Steele, J. H. and Henderson, E. W. (1995) Predation control of plankton demography. ICES J. Mar. Sci., 52, 565–573. Strang, G. (1988) Linear Algebra and its Applications. Harcourt, Brace Jovanovich, Orlando, 505 pp. Strom, S. L., Miller, C. B. and Frost, B. W. (2000) What sets lower limits to phytoplankton stocks in high-nitrate, low-chlorophyll regions of the open ocean? Mar. Ecol. Prog. Ser., 193, 19–31. Tian, R. C. C. (2006) Toward standard parameterizations in marine biological modeling. Ecol. Model., 193, 363– 386. Verity, P. G. (1991) Measurement and simulation of prey uptake by marine planktonic ciliates fed plastidic and aplastidic nanoplankton. Limnol. Oceanogr., 36, 729–750.

Murray, A. G. and Parslow, J. S. (1999) The analysis of alternative formulations in a simple model of a coastal ecosystem. Ecol. Model., 119, 149 –166.

Volterra, V. (1926) Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Mem. R. Accad. Naz. dei Lincei. Ser. VI, 2, 31– 113.

Myerscough, M. R., Darwen, M. J. and Hogarth, W. L. (1996) Stability, persistence and structural stability in a classical predatorprey model. Ecol. Model., 89, 31–42.

Wainwright, T. C., Feinberg, L. R., Hooff, R. C. et al. (2007) A comparison of two lower trophic models for the California Current system. Ecol. Model., 202, 120– 131.

Nayfeh, A. H. and Mook (1979) Nonlinear Oscillations. Wiley-Interscience, New York.

Wroblewski, J. S. (1977) A model of phytoplankton plume formation during variable Oregon upwelling. J. Mar. Res., 35, 357– 394.

1230

W. C. GENTLEMAN AND A. B. NEUHEIMER

j

FUNCTIONAL RESPONSES AND ECOSYSTEM DYNAMICS

Wroblewski, J. S. (1989) A model of the spring bloom in the NorthAtlantic and its impact on ocean optics. Limnol. Oceanogr., 34, 1563– 1571. Wroblewski, J. S. and Richman, J. G. (1987) The nonlinear response of plankton to wind mixing events—implications for survival of larval northern anchovy. J. Plankton Res., 9, 103 –123.

Wroblewski, J. S., Sarmiento, J. L. and Flierl, G. R. (1988) An ocean basin scale model of plankton dynamics in the north Atlantic 1. Solutions for the climatological oceanographic conditions in May. Global Biogeochem. Cycles, 2, 199 –218.

1231