399

Modelling copepod development: current limitations and a new realistic approach W. C. Gentleman, A. B. Neuheimer, and R. G. Campbell Gentleman, W. C., Neuheimer, A. B., and Campbell, R. G. 2008. Modelling copepod development: current limitations and a new realistic approach. – ICES Journal of Marine Science, 65: 399 – 413.

To predict the influence of environmental variability on copepod dynamics and production, models must account for the effects of temperature and food on stage-dependent time-scales. Here, data for development-time means and variance of Calanus finmarchicus are used to quantify the limitations of existing models. Weight-based individual models are sensitive to uncertain parameters, such as moulting weights, assimilation efficiency, and environmental dependencies, making them highly difficult to calibrate. The accuracy of stage-based population models using ordinary differential equations depends on model structure, with some predicted generation times being incorrect by months. Even when large numbers of age classes are used to reduce modelled variability, it is not possible to make variability consistent with the data. Accuracy of mean times for stage-based population models using difference equations requires a small time-step, which results in large numbers of age classes and modelled variability that is underestimated by orders of magnitude, unless a probabilistic moult fraction is used. We present a new stage-based individual model that avoids the limitations of other models and successfully represents C. finmarchicus mean development timing and associated variability. This approach can be adapted easily for other species, as well as dynamic environmental conditions. Keywords: copepods, development rates, food effects, IBMs, modelling, moulting, population dynamics, stage-based models, structured models, temperature effects. Received 7 July 2007; accepted 6 February 2008 W. C. Gentleman, A. B. Neuheimer: Department of Engineering Mathematics and Internetworking, Dalhousie University, 1340 Barrington Street, Halifax, Nova Scotia, Canada B3J 1Y9. R. G. Campbell: Graduate School of Oceanography, University of Rhode Island, Kingston, RI 02881, USA. Correspondence to W. C. Gentleman: tel: þ1 902 4946086; fax: þ1 902 4231801; e-mail: [email protected]

Motivation Copepods play key ecological roles in linking primary production and higher trophic levels. These roles are governed by demographic structure because copepod behaviours, physiological rates, prey, and predators change as they develop through a series of life stages. Stage-dependent time-scales, including stage duration (SDi) and time of development from egg to stage i (DTi, e.g. DTCVI, generation length), are longer and more variable when food and temperature are low (e.g. Landry, 1975; Vidal, 1980; Corkett et al., 1986; Peterson and Painting, 1990; Campbell et al., 2001). Therefore, environmental variability can alter critical population metrics, such as the average length of the non-feeding and larval phases, time between generations, time before entering dormancy, etc. Variability in development times may afford a selective advantage to individuals that can reduce or optimize these times with respect to environmental conditions (Miller and Tande, 1993; Fiksen, 2000). Thus, our ability to address many scientific questions about copepod population dynamics and production—including the effects of climate change—requires an accurate characterization of ontogenic timing. However, the mathematical models used to address these questions may not meet this requirement. Modelled development times are derived from the equations and parameters used to describe copepod dynamics, and these

# 2008

vary widely (Carlotti et al., 2000). Some models are framed in terms of populations, others individuals; some abundance, others weight; some differential equations, others difference equations, etc. This is in addition to differences regarding which processes are represented and what mathematical functions are used to describe them. For each function, parameters must be specified, and often, these values are based on laboratory experiments. Rarely are models tested for their ability to reproduce measured development times when simulating those same laboratory conditions such that their accuracy is unknown. Therefore, erroneous conclusions about factors controlling dynamics in the field are possible; modellers may attribute to other causes what is really an inherent problem in the mathematics. Here, data from a recent experiment for Calanus finmarchicus (Campbell et al., 2001) are characterized and used to evaluate quantitatively the limitations of common approaches. We also present a new model that can accurately represent laboratory data for a wide range of conditions, and is easily adaptable to dynamic field environments, as well as other species.

Experimental data Typically, measurements of copepod development times are based on experiments in which a large cohort of eggs is reared to adulthood under controlled temperature and food conditions.

International Council for the Exploration of the Sea. Published by Oxford Journals. All rights reserved. For Permissions, please email: [email protected]

400 Samples are taken at frequent intervals to assess the evolving stage structure of the population. These data provide a time-series of the percentage of the cohort that has matured to at least stage i (PCi), from which average development times and their variation are derived. The median development time of stage i (MDTi) is defined as the length of time from spawning until PCi equals 50% (Landry, 1975; Peterson and Painting, 1990; Campbell et al., 2001). Thus, MDTC1 is the average naupliar development time and MDTCVI is the average time to reach adulthood and, therefore, an estimate of generation time. The difference between the MDT of two successive stages provides a measure of the stage duration, SDi (i.e. SDi ¼ MDTiþ1 2 MDTi and MDTi ¼ SD1 þ SD2 þ . . .þ SDi21). Variation about MDTi is the result of the cumulative variation of individual spawning times and durations of stages before i, as well as sampling error. In the literature, values of MDTi have been estimated through a variety of fits to PCi (or 1 – PCi) vs. time, including least-squares linear regressions for a restricted portion of the data (e.g. 10%  PCi  90%; Peterson, 1986; Peterson and Painting, 1990; Campbell et al., 2001), sigmoidal curves (Peterson, 1986; Cook et al., 2007), and cumulative probability distributions (e.g. gamma: Klein Breteler et al., 1994; Leandro et al., 2006; Hu et al., 2007; normal: Hu et al., 2007). When different functions have been used for the same data, the resultant MDTs are nearly identical (Peterson, 1986; Hu et al., 2007), whereas the characterization of variation differs. For example, a linear fit describes development time variability (DTVi) as the inverse of the slope of the regression line for stage i (Peterson, 1986; Campbell et al., 2001), which is equivalent to using the linear fit to estimate the time interval during which the entire cohort completes the moult to stage i (Figure 1). This measure thereby neglects the development time extremes and may underestimate true variance. Although normal vs. gamma fits assume symmetrical vs. asymmetrical distributions, asymmetrical fits are not required to describe asymmetry in the data, because asymmetry arises from unequal variance in the durations of consecutive stages (Hu et al., 2007). Normal distributions arguably afford more intuitive measures, in that the median equals the mean time for an individual to reach stage i (i.e. MDTi ¼ DTi), and variation is quantified through the familiar standard deviation, sDTi. Furthermore, because of their additive properties, they imply that stage durations are also normally distributed with mean ffi stage duration SDi ¼ DTiþ1 2 DTi, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and sSDi ¼ s2DTiþ1  s2DTi , when stage durations are assumed independent. Traditional approaches to estimate probability distribution parameters require sophisticated software packages and may be prone to biases from outliers. Here, we present a simple alternative, which is akin to robust statistical approaches (Hoaglin et al., 2000) and essentially converts the “language” of the linear fit to that of a normal distribution. For a least-squares line of best fit with y-intercept a and slope b, MDTi ¼ (0.52 a)/b and DTVi ¼ 1/b. Thus, an estimate of DTi ¼ MDTi. As PCi ¼ 84% when development time ¼ MDTi þ sDTi, the regression line can be used to find t84% ¼ (0.84 2 a)/ b, which yields an estimate for sDTi ¼0.34 DTVi. The identical estimate for sDTi arises using the time at which PCi ¼ 16%, because both the regression line and normal curve are symmetrically distributed about MDTi. In other words, this formula for sDTi assumes that the regression line is a reasonable predictor of the interval over which 68% of the population moults. Estimates of sDTi are fairly insensitive to this assumption, because predicted

W. C. Gentleman et al. values change very little, even when the data are linear over smaller or larger ranges (i.e. sDTi varies by ,20% when the regression line is reasonable for a time for which anywhere from 25% to 90% of the population moults). MDT and DTV data from a recent comprehensive laboratory study (Campbell et al., 2001) were used to estimate the mean and standard deviations of development times for C. finmarchicus at three temperatures (48C, 88C, and 128C at high food) and three food levels (high  350, medium  40, and low  25 mgC l21 at 88C; Table 1). Because development time variation is cumulative, DTVi should increase with stage, but sampling and fitting error resulted in a few instances where DTVi decreased between stages (Campbell et al., 2001). Therefore, measured DTVs for stages ,CV were linearly smoothed; DTVCV were not included because their values were far greater than others (i.e. non-linearly related) and always greater than DTVCIV. Additionally, variation in spawning times was neglected in this analysis, which results in sSDEgg being an upper estimate. These data demonstrate that, at high food and 88C, DTCVI is 44 d, and roughly one-third of this generation time is spent in egg and naupliar stages (Table 1). The stage durations increase monotonically for copepodites, but not for nauplii. SDNI , SDEgg, and the stage duration of NIII, the first-feeding stage, is greater than the duration of any other naupliar stage. sDTi is small; only CVI and the non-feeding stages have coefficients of variation (C.V.; note that CV indicates the fifth copepodite stage; C.V.DTi ¼ sDTi/DTi) .10%, and all C.V.DTi are ,20%. Estimates for sDTi are relatively larger, with all C.V.SDi ,20%, except for C.V.SDCV, which is ,50%. Observations indicate that DTi varies with temperature, T, such that higher temperatures lead to shorter development times (e.g. Campbell et al., 2001). A variety of empirically derived functions have been used to describe this relationship (see Carlotti et al., 2000), with the most common being Beˇlehra´dek’s equation: DTi ¼ ai (T þ b)c. This fit implies equiproportional development (i.e. the proportion of time spent in any stage is the same at every constant temperature), resulting in the temperature dependence of stage durations also being described by Beˇlehra´dek’s equation: SDi ¼ (ai þ 1 2 ai)(Tþb)c ¼ ai0 (T þ b)c (Table 2). Using parameters for C. finmarchicus (Campbell et al., 2001), a temperature increase from 58C to 158C reduces development times and stage durations by a factor of 3 (Figure 2). Both sDTi and sSDi change very little between 88C and 128C, but exhibit increases between 48C and 88C by roughly the same proportion as DTi, resulting in little change in the C.V.s among temperatures (Table 1). One exception to this is sDTCVI, which exhibits a pronounced increase at 48C relative to 88C, which is likely the result of some fraction of the population entering dormancy (Campbell et al., 2001). Experiments have demonstrated that, when food, F, is low, development time of the feeding stages (NIII and higher) increases, even to the extent that development is arrested (e.g. Vidal, 1980; Davis, 1983, 1984; Tsuda, 1994; Campbell et al., 2001; Crain and Miller, 2001). When the food-dependent mean development times for C. finmarchicus are normalized to their temperature-dependent satiated value, the variation in food is    characterized by DTi(F) / DTi(High) ¼ (K/F)þ1 with K ¼ 17 mgC l21 (r 2 ¼ 0.97 and p ,0.01; Figure 3) for all feeding stages except NIII, because DTNIV did not exhibit any significant trend with food (Table 1). This fit is the inverse of the Michaelis– Menten ingestion rate functional dependence with a half-saturation constant of K. When food is decreased, sDTi and

401

Modelling copepod development

Figure 1. Development time of Calanus finmarchicus: circles are data from Campbell et al. (2001) for stages CIII (closed) and CIV (open) at 88C and high food (adapted from their Figure 5). Dashed line is linear regression for data 10% ,PCi ,90%, used to calculate median development time (MBTi), development time variability (DTVi), and stage duration (SDi). Curve is fit of cumulative normal distribution with DTi ¼ MDTi, sDT ¼ 0.34 DTVi.

Table 1. Mean development time (DTi) and standard deviation (sDTi) and mean stage duration (SDi) and standard deviation (sSDi) calculated from laboratory data in Campbell et al. (2001) for C. finmarchicus held under constant temperatures (48, 88, and 128C) and food conditions (high, medium, and low). Also shown is the coefficient of variation (C.V.). Temperature (88 C)/food level

i

488 C/high

888 C/high

1288 C/high

888 C/medium

888 C/low

Mean + s.d. C.V. Mean + s.d. C.V. Mean + s.d. C.V. Mean + s.d. C.V. Mean + s.d. C.V. (d) (%) (d) (%) (d) (%) (d) (%) (d) (%) DTi(d) NI 3.0 + 0.3 10 1.8 + 0.4 20 1.1 + 0.3 23 1.8 + 0.4 20 1.7 + 0.4 20 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NII 5.0 + 0.5 9 2.9 + 0.4 14 2.0 + 0.3 17 2.8 + 0.4 14 2.9 + 0.4 14 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIII 7.9 + 0.6 8 4.8 + 0.5 10 3.1 + 0.4 13 4.6 + 0.5 10 4.6 + 0.5 10 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIV 15 + 0.8 5 8.9 + 0.5 6 5.7 + 0.5 8 10 + 2.3 22 13 + 5.0 39 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NV 19 + 1.0 5 11 + 0.6 5 7.0 + 0.5 7 15 + 2.8 19 19 + 5.4 28 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NVI 23 + 1.1 5 13 + 0.7 5 8.5 + 0.6 7 18 + 3.3 18 23 + 5.9 26 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CI 27 + 1.3 5 16 + 0.7 5 10 + 0.6 6 21 + 3.7 17 27 + 6.3 23 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CII 32 + 1.5 5 18 + 0.7 4 12 + 0.7 6 26 + 4.2 16 33 + 6.7 20 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CIII 37 + 1.6 5 22 + 0.9 4 15 + 0.8 5 31 + 4.6 15 38 + 7.1 19 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . 1.8 4 26 + 0.9 4 18 + 0.9 5 37 + 5.1 14 44 + 7.5 17 CIV 45 + . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CV 56 + 2.0 3 32 + 1.0 3 22 + 0.9 4 45 + 5.6 12 58 + 7.9 14 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CVI 91 + 19 21 44 + 4.9 11 30 + 4.4 15 62 + 8.9 14 80 + 15 18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . Egg 3.0 + 0.3 10 1.8 + 0.4 20 1.1 + 0.3 23 1.8 + 0.4 20 1.8 + 0.4 20 SDi(d) . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NI 2.0 + 0.4 18 1.1 + 0.2 20 0.9 + 0.2 23 1.1 + 0.2 20 1.1 + 0.2 20 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NII 2.9 + 0.4 15 1.9 + 0.2 13 1.1 + 0.2 20 1.9 + 0.2 13 1.2 + 0.2 13 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIII 7.1 + 0.5 7 4.2 + 0.3 6 2.7 + 0.2 9 5.7 + 2.3 40 8.3 + 5.0 60 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIV 4.0 + 0.5 14 2.2 + 0.3 12 1.3 + 0.2 19 4.1 + 1.5 38 6.2 + 2.1 34 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NV 3.7 + 0.6 16 2.0 + 0.3 14 1.4 + 0.3 19 3.5 + 1.7 48 3.7 + 2.2 58 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NVI 4.2 + 0.6 15 2.5 + 0.3 12 1.8 + 0.3 16 3.3 + 1.8 54 4.4 + 2.2 51 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CI 4.9 + 0.7 14 2.8 + 0.3 11 2.1 + 0.3 14 4.9 + 1.9 39 5.4 + 2.3 43 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CII 5.7 + 0.7 13 3.5 + 0.3 9 2.4 + 0.3 13 4.9 + 2.0 41 4.9 + 2.4 48 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . 0.8 10 4.2 + 0.3 8 3.1 + 0.3 10 6.1 + 2.1 35 6.2 + 2.5 40 CIII 7.2 + . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CIV 11 + 0.8 7 6.2 + 0.3 6 4.5 + 0.3 8 8.1 + 2.2 27 15 + 2.5 17 . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CV 35 + 19 53 11 + 4.8 42 7.7 + 4.3 55 17 + 7.0 40 22 + 12 57

402

W. C. Gentleman et al.

Table 2. Empirical fit of the Beˇlehra´dek function of variation of time with temperature Time ¼ ai (T þ 9.11)22.05. Stage-dependent coefficient ai for Time ¼ DTi (originally given in Campbell et al., 2001) and ai0 for Time ¼ SDi. i ai a i0 Egg – 595 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NI 595 387 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NII 982 583 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIII 1 564 1387 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NIV 2 951 759 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NV 3 711 715 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . NVI 4 426 841 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CI 5 267 966 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CII 6 233 1 137 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CIII 7 369 1 429 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CIV 8 798 2 166 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CV 10 964 5 916 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . CVI 16 880 –

sSDi increase at rates that are proportionally greater than food effects on DTi and SDi, particularly for the late nauplii and early copepodite stages (Table 1). When compared with high food, C.V.s are increased by factors of 3 –5 at medium food, and factors of 3 –10 at low food.

Common approaches to modelling development Many mathematical models have been used to study copepod dynamics, and these models vary with respect to their units, structure, transfer functions, and parameter values. Some have been developed from a population perspective, whereas others are framed in terms of individuals, and either kind may or may not account for variation among individuals. An excellent review is provided by Carlotti et al. (2000). Here, some common approaches to modelling development are introduced, and their ability to represent the laboratory data in Table 1 is examined.

Figure 2. Variation of modelled stage durations (SDi) with temperature relative to the value at 88C. Solid line is empirical fit to data using the Beˇlehra´dek equation (Table 2). Dotted lines are for inverse Q10 relationship for Q10 values indicated. Note that curves are the same for all i.

Figure 3. Variation of modelled stage duration (SDi) with food relative to the value at high food condition, resulting from different ingestion functions in weight-based models. Solid line is empirical fit using the inverse Michaelis – Menten (MM) equation (in the section Experimental data) with Ksat ¼ 17 mgC l21. Dashed line is MM with Ksat ¼ 70 mgC l21. Dotted line is modelled functional response used in Carlotti and Radach (1996), CR96. Note that curves are the same for all feeding stages (i ¼ NIV 2 CV). For these data–model comparisons, modelled values were either generated explicitly (e.g. mean development times, in the section Weight-based models of individuals, or analytical functions for PCi, in the section Ordinary differential equations), or estimated using the same approach as was used for the data (i.e. compute modelled MDTs and DTVs, and use those to estimate means and variation for all 13 stages, as explained in the section Experimental data).

Weight-based models of individuals Weight-based models of individuals describe the time-dependent change of an individual copepod’s weight (W ), measured in terms of a particular nutrient (e.g. carbon). They are typically framed as ordinary differential equations (ODEs) for dW/dt, the net growth rate, equal to the sum of ingestion, egestion, respiration, excretion, reproduction, etc. (e.g. Steele and Frost, 1977; Steele and Mullin, 1977; Batchelder and Miller, 1989; Carlotti and Sciandra, 1989; Carlotti and Nival, 1992; Batchelder and Williams, 1995; Carlotti and Radach, 1996; Heath et al., 1997; Batchelder et al., 2002). Each physiological rate generally varies with W (i.e. is allometric) and may vary with environmental conditions, such as food and temperature. These models are essentially the individual’s analogue to the dZ/dt equation for zooplankton biomass in the common nutrient –phytoplankton –zooplankton (NPZ) ecosystem models. Weight-based models of individuals differ in terms of which processes are considered, and which functional dependencies and associated parameter values are used (see review in Carlotti et al., 2000). Representing development through the life stages in weight-based models requires a secondary set of assumptions that relate weight with stage. Typically, this is achieved by defining a set of moulting weights, so that, when copepods in stage i reach a certain size, Wi,moult, they moult to the next stage (e.g. Batchelder and Miller, 1989; Carlotti and Radach, 1996; McLaren, 1997). Thus, the ability to model development times

403

Modelling copepod development accurately depends on at least two factors: (i) characterization of Wi,moult, and (ii) parameterizations of the growth model.

Effect of specified moulting weights on mean development times The concept of “moulting weights” is supported by laboratory data that indicate no overlap of copepodite structural weights among life stages at a specific temperature (McLaren, 1986; Carlotti et al., 1993). However, structural weight transition values are not constant among studies (Table 3) and are unknown for some species, leading modellers to use a variety of approaches for designating Wi,moult. Some assume moulting weights increase by a constant proportion between stages (e.g. McLaren, 1997). Others use experimentally measured mean dry weights, which can vary from study to study even when conditions are identical (e.g. Carlotti et al., 1993 vs. Campbell et al., 2001). Still others have used weights derived from field measurements of other species (e.g. Heath et al., 1997). A simple calculation can be used to illustrate the sensitivity of modelled development times to the choice of Wi,moult. Measured growth rates often assume an exponential growth rate within a stage (Hirst et al., 2005; Kimmerer et al., 2007), implying that, when stage-dependent growth rates, gi, are known, the time spent growing from Wi21,moult to Wi,moult is SDi ¼

  1 Wi;moult : ln gi Wi1;moult

Figure 4. Modelled stage durations relative to those measured at 88C in Campbell et al. (2001) for different sets of copepodite moulting weights. (a) – (e) correspond to sets of copepodite weights in Table 3. Tot is total time in copepodite stages CI – CV. For all cases, stage-specific growth rates were equal to the average of between-stage carbon growth rates [e.g. gCI ¼ (g*NVI2CI þ g*CI2CII)/2, where g* were taken from Campbell et al. (2001), and WNVI,moult was 1.1 mgC as in Carlotti and Radach (1996)]. is illustrated with a simple growth model:

ð1Þ

Thus, modelled stage duration is set by the ratio of the weights at which a copepod moults out of a stage vs. into it. Candidate values for Wi,moult vary widely in the literature (Table 3), and these different sets result in modelled stage durations that range by 1 –20 standard deviations about the mean (Figure 4; Table 1). Total copepodite development time (i.e. time from new CI to new CVI) ranges from 23 to 36 d, which is a 25% spread around the measured value. No one type of data appears to be the best predictor: copepodite development times are underestimated using structural weights, and both underestimated and overestimated using mean field and laboratory dry weights (Figure 4).

Effect of growth model parameterizations on mean development times For any given set of moulting weights, development times can vary owing to differences in the growth model parameterizations. This Table 3. Candidate sets of moulting weights of C. finmarchicus copepodite stages from the literature. All weights measured in mgC. (a) upper structural weights, Carlotti and Radach (1996); (b) upper structural weights, McClaren (1986); (c) mean field dry weights, Heath et al. (1997); (d) mean field dry weights, Campbell et al. (2001); (e) mean laboratory dry weights, Campbell et al. ( 2001). Stage (a) (b) (c) (d) (e) CI 2.5 5 2.0 1.4 3.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . CII 7.0 12 4.2 4.3 7.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . CIII 15 25 10 12 19 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . CIV 40 60 19 40 48 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. . . . . CV 90 150 48 136 280

dWi ¼ gi Wi1 : dt

ð2Þ

The term on the right is used for ingestion and assimilation in most weight-based models and is the form of net growth in many models (e.g. Batchelder and Miller, 1989; Batchelder and Williams, 1995; Batchelder et al., 2002), or growth when basal metabolic losses are negligible (e.g. Carlotti and Sciandra, 1989; Carlotti and Radach, 1996). The parameter 1 usually ranges between 0.6 and 1.0, with 0.7 –0.75 being the common value for copepods (van den Bosch and Gabriel, 1994; Carlotti et al., 2000). gi is a variable that is derived from a number of factors including assimilation efficiency, temperature, and food-limitation effects. For constant environmental conditions, gi has a fixed value, so that Equation (2) yields SDi ¼

1 11 Þ; ðW 11  Wi1;moult gi ð1  1Þ i;moult

when

1 = 1:

ð3Þ

There is a non-linear relationship between 1 and the stage durations that results in any uncertainty in 1 being amplified for SDi. For example, with WCIV,moult ¼50 mgC and WCV,moult ¼100 mgC, the simple choice of 1 ¼ 0.7 vs. 1 ¼0.75—a difference of 7%— results in modelled stage durations for CV that differ by 20%. The inverse relationship between gi and SDi results in any change in gi introducing a change of a reciprocal amount in SDi (e.g. halving gi doubles SDi). For example, assimilation efficiency usually ranges between 0.6 and 0.8 (Carlotti et al., 2000), which results in the corresponding modelled SDi ranging by .40%. Of the weight-based models that include temperature effects, the same relationship is usually applied to all rates and stages (e.g. Batchelder and Williams, 1995; Carlotti and Radach, 1996). Temperature dependence is often described using a Q10 value,

404 which is the relative amount a physiological rate, r, increases if the temperature rises by 108C above a reference temperature, Tref (i.e. Q10 ¼ r(Tref þ 10)/r(Tref ). It is formulated as r(T ) ¼ ref )/10 r(Tref )Q (T2T , which is mathematically identical with the 10 exponential formulation r(T ) ¼ r(Tref )eq(T2Tref ), where q ¼ ln (Q10)/10. Thus, a Q10 effect on the growth rate would alter gi, and the corresponding change in SDi would be SDi,model (T ) ¼ SDi,model (Tref )e 2q(T2Tref ). There is no significant difference between the inverse Q10 dependence of stage durations and the Beˇlehra´dek equations presented in the section Experimental data when the parameterizations are chosen to match observed variation in the development times (Figure 2). However, the Q10 used in the literature is based on growth (2–2.5: Batchelder and Williams, 1995; Carlotti and Radach, 1996; Heath et al., 1997; Carlotti et al., 2000), and is typically lower than the empirical value for development (3.6: Campbell et al., 2001), resulting in development at cold temperatures being poorly characterized. For example, at 48C, a Q10 of 2 underestimates stage durations of C. finmarchicus by 30%, which would underestimate the generation time by almost 3 weeks. In virtually all weight-based models, ingestion rates—and other rates tied to ingestion, e.g. egestion—vary with food, F, so that growth increases as F increases. This relationship is usually specified by the Michaelis –Menten equation [e.g. Ii ¼ (F/(Ksat,i þ F)) Imax,i . . ., where Ksat is the half-saturation constant, and Imax the maximum ingestion rate], which is equivalent to a Holling Type 2 functional response (Gentleman et al., 2003), although other formulations are also used. However, modelled dependencies often do not represent observations. For example, the Michaelis– Menten equation approaches satiation very slowly, such that ingestion rates are only 67% of their maximum at F ¼ 2Ksat, and 83% at F ¼ 5Ksat. Thus, for a typical Ksat ¼ 70 mgC l21 (e.g. Batchelder and Williams, 1995), copepod growth is modelled as moderately food-limited, even at bloom-level phytoplankton densities. Other formulations also characterize observed low food responses poorly, such as no ingestion for F , 50 mgC m23 (Carlotti and Radach, 1996), whereas data indicate positive growth rates for food levels half that amount (Campbell et al., 2001). Because the modelled effect of food limitation on growth is exaggerated (i.e. gi is decreased by too much), the corresponding lengthening of development times is also exaggerated. For example, at 100 mgC l21, development times should increase by 17% relative to their satiated values, whereas the functional response used by Carlotti and Radach (1996) predicts an increase of 200% (Figure 3). This model also predicts arrested development at the medium and low food conditions, which was not observed. When Michaelis –Menten with Ksat ¼ 70 mgC l21 is used, modelled development times are more than double their food-satiated value at F ¼ 50 mgC l21, as opposed to the actual increase of 34% (Figure 3). This parameterization also results in a poor representation of development at high food levels, overestimating development time by 23% at satiating conditions (Figure 3).

A note on variability of timing in bioenergetics models Most weight-based models do not include any stochasticity, i.e. their parameter values are identical for all individuals. Thus, a modelled population of individuals subject to identical conditions will form a tight cohort that grows at the same rate and will moult in synchrony when they reach their moulting weights. In this case, the modelled DTV will be the model time-step, which is typically ,1 d, and the models will far underestimate observed variability.

W. C. Gentleman et al. Some models have been designed to include explicitly variable growth and initial weights (McLaren, 1997) or distributed moulting weights (Carlotti and Radach, 1996). However, there is no straightforward connection between the specific distributions used and those observed for development, making calibration difficult.

Stage-based models of populations Stage-based models of populations are designed to simulate the demographic structure by dividing the total abundance into categories that differ in size, behaviour, prey, predators, etc. The term “stage” is a misnomer in that the modelled categories are not always the same as the 13 life stages of copepods. For example, some models aggregate life stages, considering early vs. late nauplii and copepodites (e.g. Lewis et al., 1994; Plaga´nyi et al., 1999; Soussi and Ban, 2001; Hu et al., in press). Others use more than one model stage to simulate different phenologies of a particular life stage, such as having two CV stages to account for diapausing vs. active behaviour, or three adult stages to account for males vs. reproductively immature females vs. reproducing females, or for tracking life stages from different generations (e.g. Lynch et al., 1998; Plaga´nyi et al., 1999). Additionally, most models subdivide the stages into age classes as discussed below. Although the term “group-based” population models is arguably more appropriate, “stage-based” is now established in the vernacular and will be used here. The general idea of these models is to solve for the temporal variation of Ni,j, the number of animals in stage i and age class j, where there are a total of n distinct stages (i ¼ 1 to n) with each stage made up of a total of mi age classes ( j ¼ 1 to mi). Constraining our discussion to the terms related to development (i.e. ignoring egg production, mortality, etc.), changes in Ni,j are caused by (i) the ageing of animals from age class j to j þ 1, and (ii) the moulting of animals from stage i into the first age class of i þ 1 (Figure 5). For a given stage, the rate of ageing relates to the age-class duration, ti, and the rate of moulting is related to stage duration, SDi. Although some models calculate SDi using weight-based models of individuals (e.g. Carlotti and Sciandra, 1989; Carlotti and Nival, 1992; Carlotti and Radach, 1996), the traditional stage-based population models discussed here do not include variables for weight, and instead specify SDi directly (i.e. growth and development are decoupled). Previous authors have noted that use of a single age class (i.e. mi ¼ 1 for all i) vs. multiple age classes (mi . 1) affects the ability to represent development in stage-based population models accurately (e.g. Sciandra, 1986; Soussi and Ban, 2001). However, as we demonstrate here, other aspects of these model structures are equally, if not more, important. Modelled development times are influenced by the mathematical form of the equations, with the two most common being ODEs, wherein instantaneous rates of change are specified directly, and difference equations, which describe the overall change in numbers over a finite period of time. Modelled times also depend on the number of age classes considered and whether or not moulting occurs from a single age class or in a distributed manner (Figure 5a vs. b).

Ordinary differential equations Traditional stage-based models of development using ODEs describe dNij/dt, the instantaneous rate of change of Ni,j as:

405

Modelling copepod development Ageing: dNdt1;1 ¼ r1 N1;1 dNi;j ¼ ri Ni;j1  ri Ni;j dt

i ¼ 1 to n;

j ¼ 2 to mi

ð4Þ

Moulting (and ageing): dNiþ1;1 ¼ ri Ni;mi  riþ1 Niþ1;1 dt

i ¼ 1 to n  1;

where ri is the rate of maturation (ageing or moulting) of that stage and age class [T 21], often called the development rate. Note that moulting from stage i to stage i þ 1 occurs only after the copepods have matured through the entire series of mi age classes (Figure 5a). The age-class duration ti ¼ SDi/mi, meaning ti are (i) generally not the same among stages because SDi do not differ by integer amounts, and (ii) not constant in dynamic environments for which SDi varies in time because mi is a fixed aspect of the model. The coupled differential equations are subsequently integrated—typically with a numerical method (e.g. Euler or Runge–Kutta) with model time-step, Dt ¼ minutes to hours—to solve the temporal evolution of Ni,j. The ability of ODEs to represent development time means and variability accurately lies in the choice of values for ri, n, and mi, which are not independent choices. A range of parameterizations has been used to relate ri to the age-class duration ti, including ri ¼ 2/ti (Wroblewski, 1980, 1982), ri ¼ 1/ti (Sciandra, 1986; Lewis et al., 1994; Lynch et al., 1998), and ri ¼ (ln 2)/ti (Zakardjian et al., 2003). To illustrate the effect of these different parameterizations, first consider a

simple scenario where the population is partitioned into two model stages each, with only one age class (i.e. n ¼ 2, m1 ¼ m2 ¼ 1, and no j index is required). Thus, N1, the abundance in the first model stage, is the total population in stages egg to i, and t1 ¼ SD1 ¼ DTi. When simulating the maturation of a single cohort under constant laboratory conditions, t1 and therefore r1 are constant and Equation (4) yields PCi;model ðtÞ ¼

  N2 ðtÞ ¼ 1  er1 t : ðN1 ðtÞ þ N2 ðtÞÞ

ð5Þ

The PCi,model curves resulting from the different parameterizations of r1 are all similarly shaped, with larger r1 exhibiting larger initial slopes (Figure 6). However, this shape is very different from the data. The model neither exhibits a time delay before development nor a small interval of ages over which the organisms moult. Instead, development begins immediately and occurs over a range of ages well beyond what is observed, resulting in modelled coefficients of variation (C.V.DTi,model) 100%—more than an order of magnitude larger than that observed in the data (Table 1). Note that, although the choice of r1 ¼ ln 2/t1 results in DTi,model ¼ DTi, this choice also results in the worst modelled variability (Figure 6). The two other choices of r1 have somewhat improved, sDTi,model, but their DTi,model underestimate DTi by 65% and 30%, respectively (Figure 6). Therefore, it is not possible simply to adjust the parameter values to match the model and data; the underlying mathematics need to be altered. The problem outlined above is sometimes referred to as numerical dispersion and is well known in the literature. It arises because use of a constant development rate for a stage is

Figure 5. An illustration of the structure of stage-based models of populations where Ni,j is the number of animals in stage i and age class j, which changes over time owing to ageing and/or moulting. (a) Moulting occurs only from a single age class; and (b) moulting is distributed such that multiple age classes contribute to stage i þ 1.

406

W. C. Gentleman et al.

Figure 6. Effect of rate parameterizations on stage-based population models using ODEs. Plots of cumulative percentage of an initial cohort of newly spawned eggs that has reached at least stage i (PCi) vs. time since spawning, which result from a model using ODEs for two stages with a single age class [Equation (5)]. Curves for different rate parameterizations are shown, and t1 ¼ DTi. Also shown is characteristic PCi, data (solid line). effectively the same as assuming a uniform distribution of ages within that stage, which is clearly not the case when simulating a maturing cohort. A simple fix is achieved by subdividing the stage (i.e. creating mi . 1 age classes), which allows for varying age-within-stage distributions. The rate of ageing is modelled as the reciprocal of the age-class duration, t1, so that the dynamics within each age class will be like those in Figure 6, with only 68% of the animals maturing through the age class within t1. However, use of the subdivisions results in the overall moulting rate for the stage to be changed from a simple exponential to a sum of exponentials, thereby allowing for a time-varying moulting rate. For example, if we subdivide stage 1 of our two-stage model into m1 age classes so that now t1 ¼ SDi/m1 ¼ DTi/m1,

PCi;model ðtÞ ¼ 1 

m1 X j¼1

er1 t

ðr1 tÞ j1 : ð j  1Þ!

Figure 7. Effect of subdivisions on stage-based population models using ODEs. Curves for two-stage model where group 1 has m1 age classes with duration ti (ri ¼ DTi/m1, ri ¼ 1/ti), with m1 varying as shown. underestimated by 50% or overestimated by 40% (Figure 8, Columns b and c). Age classes are necessary to reduce modelled variability to the observed levels, and large numbers may be required. Furthermore, it may not be possible for modelled variability to be consistent with the data, because C.V.DTi,model decreases as the number of younger ages and/or earlier stages increase (Figure 7), whereas C.V.DTCVI is typically larger than most other stages (Table 1). For example, for a 13-stage model of C. finmarchicus at 88C and High food with ti , 0.2 d for every stage (i.e. 225 variables), sDTi,model are still overestimated by factors of 2 –2.5 for all copepodite stages except CVI, which is underestimated by a factor of 2 (Figure 8, Column d). In addition, there is no obvious way to

ð6Þ

When r1 ¼ 1/t1, DTi,model rapidly converges on DTi as m1 increases, so that the error is ,5% for ten age classes (Figure 7). Increasing m1 also reduces the modelled variability. However, even with 100 age classes, C.V.DTi,model ¼ 10%, which is still larger than the data for most stages (Table 1). It must be emphasized that, as far as the mathematics are concerned, increasing the demographic resolution through subdivision of a stage is akin to increasing the number of stages younger than i. That is, using n . 2 has a similar effect on DTi,model and C.V.DTi,model, as is shown for mi .1 in Figure 7. A model with 13 stages, each with a single age class, and ri ¼ 1/ti, results in DTi,model that are accurate to ,3%, with C.V.DTi,model ¼ 30 –40% for copepodite stages (Figure 8, Column a). This fact appears to be poorly understood in the literature, because stagebased models are accused of being unable to represent developmental delays (Soussi and Ban, 2001) and are incorrectly parameterized for ri, such that modelled generation times are

Figure 8. Modelled development times (DTCVI) and standard deviations (sDTCVI) for stage-based models using ODEs compared with data for stage-based population dynamics using differential equations (in the section Ordinary differential equations). Shown are results for generation time (i ¼ CVI) for a 13-stage model with (a) a single age class (ti ¼ SDi), and ri ¼ 1/ti; (b) a single age class and ri ¼ 2/ti (Wroblewski, 1980); and (c) a single age class and ri ¼ ln 2/ ti (Zackardjian et al., 2003); and (d) multiple age classes (ti , 0.2 d for all i). Data are for 88C and high food conditions as in Table 1.

407

Modelling copepod development relate development time variability to environmental conditions when these are dynamic, because mi is a fixed aspect of the model structure.

Difference equations Traditional stage-based population models, formulated with difference equations, use values at time t to characterize the integrated effects of ageing and moulting on Ni,j at a future time t þ Dt, where typically Dt ¼ 1 d (Davis, 1984; Miller and Tande, 1993; Plaga´nyi et al., 1999; Hu et al., 2007). Ageing is a discrete process, whereby the age-class duration ti ¼ Dt, and 100% of the animals in stage i and age class j that do not moult during the time-step are transferred to the j þ 1 age class. In their most general form, moulting occurs from a distributed range of age classes (Figure 5b), so that Ageing:   Ni;jþ1 ðt þ DtÞ ¼ 1  fi;j ðtÞ Ni;j ðtÞ

i ¼ 1 to n; j ¼ 1 to mi1

ð7Þ

Moulting: Niþ1;1 ðt þ DtÞ ¼

mi X

fi;j ðtÞNi;j ðtÞ

i ¼ 1 to n  1;

j¼1

where fi,j (t) is the fraction (0  fi,j (t)  1) of the population in stage i and age class j that moult to stage i þ 1 between times t and t þ Dt. These coupled equations can be written in matrix form, so that this is one type of matrix population model (e.g. Caswell, 1989). The difference equations are applied repeatedly, stepping forward in time, to generate a series of Ni,j at discrete times. Because the numerical solution of differential equations effectively converts them to difference equations (e.g. forward Euler: dNi,j(t)/dt  Ni,j (t þ Dt) 2 Ni,j (t)/Dt), solution of the ODE models (in the section Ordinary differential equations) with a single age class (i.e. mi ¼ 1 for all i, no j subscript required) also results in Equation (7) with fi ¼ riDt and Dt ¼ minutes. Thus, the numerical dispersion issues discussed previously apply for difference equation models, indicating that multiple age classes are needed to represent development-time variability accurately. Note that ODEs and difference equations are fundamentally different when either (i) Dt is sufficiently large that the derivative is poorly approximated (e.g. for a 13-stage, single age-class model, errors in DTi,model are doubled for Dt ¼ 1 vs. 0.01 d); or (ii) multiple age classes are used, because traditional ODEs do not represent ageing discretely or have the same age-class duration for all stages. Thus, difference equation models present their own unique challenges to representing development times and variability accurately, and their ability to do so relates to the formulation of fi,j and the size of Dt. The most intuitive formula for fi,j is a switch where the copepods spend a period equal to at least SDi in the stage before moulting (Figure 5a), i.e.  fi;j ¼

Although this approach seems rational when considering constant laboratory conditions, it is not appropriate when modelling a dynamic environment wherein SDi varies in time. As a simple illustration, consider a situation where environmental conditions were constant, the corresponding SDi ¼ 10 d, and a copepod was in that stage for 5 d. Now suppose conditions changed so that the stage duration for the new conditions was 5 d. One would not expect the copepod to moult suddenly because its age now corresponded to the new stage duration; it had only been in that stage for a period corresponding to 50% of its development. This leads to the concept of the moult cycle fraction (MCF; Miller and Tande, 1993), which measures proportional development such that MCFi ¼ 0 when a copepod first moults into stage i, and it increases over time until MCFi ¼ 1, at which point the copepod moults to stage i þ 1. The amount by which MCFi increases between t and t þ Dt depends on the development rate during that time, which is approximated by the reciprocal of the stage duration, i.e. MCFi,jþ1 ¼ MCFi,j þ Dt/SDi. Thus, MCF is effectively the discrete analogue of the variable age-class duration inherent in ODE models (in the section Ordinary differential equations). With the MCF approach, copepods will moult when their age equals the stage duration in a constant environment, but development time is adjusted appropriately when conditions are varying. Using the example in the previous paragraph, MCFi,j ¼ 0.5 just before the change. On the next day, MCFi,j+1 ¼ 0.5 þ 1/5 ¼ 0.7. If the environment remained constant at the new conditions, it would take an additional 1.5 d to moult. In other words, the total time spent in the stage would be 7.5 d, an average of the stage durations at the two conditions. Thus, for dynamic environments, the correct formulation for a switch is  fi;j ¼

MCFi;j , 1 MCFi;j  1:

0; 1;

Formulations such as Equations (8) and (9) result in a cohort moulting in synchronicity. That is, at the first time-step when MCFi,j  1, 100% of a particular cohort will moult, resulting in SDi,model ¼ age(MCFi,j  1) 2 Dt/2. Because the age increments and model time-step are identical, the chosen time-step limits the ability to resolve the stage duration (Miller and Tande, 1993). For example, for C. finmarchicus at 88C and high food, SDNI ¼ 1.1 d (Table 1); however, when Dt ¼ 1 d, copepods will not moult until their age ¼ 2 d, and SDNI,model will have an error of 0.4 d. Similar age discretization error will be incurred at each moult, with the relative importance being greater for stages with SDi closer to Dt (Hu et al., in press). The associated accumulated error in DTi will depend on the signs and sizes of discrete age error incurred in earlier stages, which can be significant. For a 13-stage model of C. finmarchicus at 88C and high food with Dt ¼ 1 d, errors DTi,model in range from 15% to 26% among stages (Figure 9, Column a). It is possible to correct for this discrete age problem by keeping track of the fraction of the time-step spent in the new stage after moult, and assigning the corresponding amount to MCFiþ1,1, i.e. 

0; 1;

age ðjÞ , SDi age ðjÞ  SDi :

ð8Þ

ð9Þ

MCFiþ1;1 ¼

 MCFi;j  1 SDi : SDiþ1

ð10Þ

408

W. C. Gentleman et al. not moulted by age class j. As the difference between the cumulative and conditional probabilities increases with variance, errors associated with use of fi,j ¼ Pi,j depend on the relative standard deviations. For a 13-stage model of C. finmarchicus using distributions based on the data in Table 1 at 88C and high food and Dt ¼ 0.1 d, DTi,model is accurate to within 5% for all stages except for egg (9%) and CVI (20%; Figure 9, Column c). Hu et al. (2007) proposed a “corrected” probability formula: fi;j ¼

Figure 9. Modelled development times (DTCVI) and standard deviation(sDTCVI) compared with data for stage-based population dynamics using difference equations (in the section Difference equations). Shown are results for generation time (i ¼ CVI) for a 13-stage model with multiple age classes. (a) f ¼ switch with MCF ¼ 0 at moult as in Equation (8), Dt ¼ 1 d; (b) as in (a) but with Dt ¼ 0.1 d; (c) f ¼ probabilistic model with cumulative normal (Davis, 1984), Dt ¼ 0.1 d; (d) f ¼ probabilistic model as in Equation (10) with Dt ¼ 0.1 d. Data are for 88C and high food conditions as in Table 1. Equation (10) is effectively what Hind et al. (2000) do in their formula for cumulative development [their Equation (10)], although they do not account for the ratio of stage durations, but the error introduced by this oversight is quite small even when Dt ¼ 1 d. Application of Equation (10) is limited because (i) it cannot be extended to situations in which more than one group contributes new moults to stage i þ 1 during a time-step (e.g. moults occur from multiple age classes owing to individual variability); and (ii) even with Dt ¼ 1 d, it suffers from the problem that the modelled variability is linked to the time-step and is still too small (i.e. sDTi,model ¼ 0.34 d), especially for food-limited conditions. A more practical and common way to deal with the discrete age problem is to reduce the time-step so that the discrete age error is negligible and animals effectively have MCF ¼ 0 at moult. For C. finmarchicus at high food and 88C, a time-step of 0.1 d is needed to have reasonably accurate DTi,model (errors ,5% for all 13 stages; Figure 9, Column b). However, this results in sDTi,model ¼ 0.034 d, which is orders of magnitude too small (Figure 9, Column b). The natural way to correct modelled variability is to make development distributed, so that more than one age class can moult to the next stage (Figure 5b). Mathematically, this is done by using an fi,j that is probabilistic, and in fact, many modellers call fi,j the “probability of moulting” for stage i and age class j (Davis, 1984; Carlotti et al., 2000; Hu et al., 2007). However, use of that term is somewhat misleading because it does not indicate which kind of probability, causing some modellers to use erroneous formulations that introduce error into the modelled mean times, not just the variance. For example, Davis (1984) used fi,j ¼ Pi,j, the cumulative probability for age j relative to stage i’s mean stage duration and standard deviation. However, Pi,j is the probability that the original population moults by age class j, whereas fi,j relates to the probability of moulting for animals remaining at age class j (Hu et al., 2007). In other words, fi,j is the conditional probability that individuals will moult by age class j þ 1, given that they have

pi:j ; 1  Pi:j

ð11Þ

where pi,j is the probability density function corresponding to the cumulative probability Pi,j. Although Equation (11) is the correct form for the instantaneous moulting rate, it is not correct for use as fi,j in the difference equation [Equation (7)], because that fraction is an integrated quantity. Equation (11) is erroneous by a factor Dt and, when employed in our model, results in fi,j . 1 so that more modelled animals moult than exist in the population. The correct formula is Pi;jþ1  Pi;j  ; fi;j ¼  1  Pi;j

ð12Þ

which results in accurate DTi,model and si,model with ,2% error for all stages (Figure 9, Column d) when Dt ¼ 0.1 d. The error in the formula proposed by Hu et al. [2007; i.e. Equation (11)] was not apparent because their choices of (i) Dt ¼ 1 d, and (ii) approximating pi,j using piecewise linear functions resulted in their particular implementation of Equation (11) being exactly the same as Equation (12). Two additional issues arise with probabilistic moulting. The first relates to the choice of mi, because probabilistic moulting allows some copepods to have ages .SDi (Figure 5b). When the total time allotted to a stage ¼ SDi þ 3sSDi, the number of variables increases by .50%, compared with those required by the switch formulation when Dt ¼ 0.1 d. Additionally, there will be some copepods that do not moult before reaching the last age class, and a modelling decision has to be made about what to do with these non-moulters (e.g. have them die). The second issue, the appropriate way to alter Equation (12) for a dynamic environment, is not obvious, because (i) the formula requires cumulative probabilities, which are based on constant conditions; and (ii) it is nonsensical to use a probability distribution based on MCF because, by definition, animals moult when MCF ¼ 1.

A note on hybrid model structures Other stage-based population models have been formulated using a hybrid approach, wherein ageing is described using difference equations and moulting is described by ODEs. In their original incarnation, hybrid models characterized instantaneous moulting rates, using modelled weights (e.g. Carlotti and Sciandra, 1989; Carlotti and Nival, 1992; Carlotti and Radach, 1996), but later versions decoupled growth and development by replacing the weightbased rates with age-based probabilities (Soussi et al., 1997; Soussi and Ban, 2001), although the incorrect probability formulation was used (i.e. cumulative vs. conditional; Hu et al., 2007). This results in distributed moulting (Figure 5b), with ODEs solved over smaller time-steps than the age-class duration (e.g. Dt ¼ 1 h vs. t ¼ 6 h). It is not immediately apparent whether there are advantages to this hybrid approach over traditional difference

409

Modelling copepod development equation models that would justify the increased computational cost. Furthermore, this method suffers from the same limitations as difference equations in that large numbers of variables are required to represent development accurately, and the appropriate way to alter cumulative probabilities for a dynamic environment is uncertain.

at the same deviation from the mean for every stage: Kimmerer et al., 2007; or updated at every time-step: Hu et al., 2007) cannot reproduce the laboratory data standard deviations even with small Dt.

An alternative approach: stage-based models of individuals

Here we present a simple method of estimating means and standard deviations of stage-dependent time-scales from standard laboratory measurements. These data were provided for C. finmarchicus for a range of experimental conditions along with empirical fits for the environmental dependencies (Campbell et al., 2001). Both temperature and food have pronounced effects such that development times increase by roughly 50% for a 48C decrease in temperature or a significant reduction in food (Tables 1 and 2: 88C/high vs. 48C/high and 88C/low). Thus, for typical ranges of field conditions during the copepods’ active period (e.g. 2 –128C and 0 –350 mgC l21 for winter/spring on Georges Bank), mean development times will vary by factors of 2–3 (Table 1) and potentially much more owing to the multiplicative effect of cold temperatures before the spring phytoplankton bloom. Food effects may also mitigate temperature-related changes and vice versa, if one factor increases while another decreases. Food exerts a stronger effect on development-time variability than temperature, such that C.V.DT at low food are 4–8 times those at cold temperatures. Considering that 95% of the population moults within the window DTi þ 2sDT, there will be generational overlap even at moderate food levels. Furthermore, because the spread in development times is increased at low food, populations may be buffered from effects of food variability (e.g. shifts in timing of spring phytoplankton bloom), with the buffering extent depending on concurrent temperature variability. Without modelling, it is difficult to know a priori whether temperature or food is more important in controlling metrics such as naupliar development and generation times, especially because the frequency of variability of these two factors is generally different. Temperature usually changes by a few degrees Celsius each month, whereas food variability is often approximately days to weeks. Even daily, variability of both temperature and food can be quite large, as when copepods migrate vertically through stratified layers. Therefore, studying the influence of climate on copepod population dynamics and production requires models to have accurate characterization of the effect of both temperature and food on stage-dependent time-scales in addition to their associated variability. As demonstrated here, many existing models are limited in their ability to predict development time means and variability even under constant environmental conditions, let alone dynamic field environments.

There is another rarely used approach to modelling copepod development: stage-based models of individuals (IBM; Miller et al., 1998; Crain and Miller, 2001; Hu et al., 2007; Kimmerer et al., 2007). Like the weight-based models, these explicitly consider each individual. However, instead of using equations for growth, they are formulated to represent the development through stages directly, like the stage-based models. An individual is assigned an initial stage and age (or MCF) and, as time progresses, the age (or MCF) increments until age equals the stage duration (or MCF ¼ 1), whereupon the individual moults to the next stage. Previous stage-based models of individual copepods were of limited application, in that (i) they were not framed for use in simulations where mean development times and individual variability can change with environmental conditions (e.g. Hu et al., 2007; Kimmerer et al., 2007); (ii) they neglected variability in that the same stage durations were used for each individual (e.g. Miller et al., 1998); and/or (iii) they either neglected food limitation effects or else approximated them by arresting development for a random fraction of the population (Miller et al., 1998; Crain and Miller, 2001). To our knowledge, no model of this type has attempted to link MCF to food effects using empirical data. Our new model builds on previous stage-based models of individuals by including laboratory data for both temperature and food effects, as well as individual variability. An initial population is established, which entails specifying the total number of individuals, and each individual’s stage i, MCFi, and Fitnessi, where Fitnessi is a uniformly distributed random variable. The corresponding stage duration of the individual, SDi, is found by using the inverse of the cumulative distribution function (e.g. norminv function in Matlab’s statistics toolbox; Mathworks, 2005) for 1-Fitnessi, SDi, and sSDi, so that higher Fitness corresponds to shorter stage durations (Figure 10). At each time-step, every individual’s MCFi increments by Dt/SDi as in the case for difference equations, but note that the increment value will vary among individuals, because of the individual-dependent stage duration. Thus, even a tight initial cohort will soon exhibit a spread in stages and age-within-stage, and will not moult synchronously. Individuals with MCFi . 1 are identified, and their stage is advanced by 1. They are assigned a new Fitnessi, which is used to calculate a new SDi. Their MCFi is then set equal to the fraction of Dt that they would have spent in the new stage divided by the new stage duration, as in Equation (10), to avoid age discretization error. Even with large Dt (e.g. 1 d), the resulting model has DTi,model accurate within ,3% for all stages, and sDTi,model ranging between 50% and 130% of the observed value for all stages of C. finmarchicus and every environmental condition (Figure 11). This result is striking given that all previous data– model comparisons were restricted to only the 88C and high food condition. Furthermore, use of a smaller Dt results in all sDTi,model accurate to within 7%, indicating that the error for Dt ¼ 1 d has to do with temporal sampling, not modelling approach. Models using alternative schemes for assigning individual variability (e.g. fixed

Discussion

Limitations of weight-based models Development times in weight-based models are sensitive to the choice of moulting weights, which are used to relate weight to stage (i.e. growth to development). Different values from the literature result in modelled stage durations that range by factors of 2 –5 among sets, and these ranges apply for all environmental conditions, because moulting weights are specified as constants (i.e. 2 weeks for CV at 88C and high food, and longer at lower temperature or food conditions). Furthermore, no one type of data (i.e. structural, field, or laboratory) appears to be the best predictor, because all types underestimated the duration of some

410

W. C. Gentleman et al.

Figure 10. The effect of incorporating individual variation in fitness by randomly sampling a uniform distribution. The individual’s stage duration is determined by using the inverse cumulative distribution function (CDF), e.g. normal CDF with higher fitness corresponds to shorter stage durations. stages and overestimated others. To improve accuracy, more research into moulting weights is needed. This includes the need to characterize environmental dependencies (e.g. copepod weight is reduced when food is low or temperature is high; Miller et al., 1977; Carlotti et al., 1993) and to determine how best to represent the critical weight concept for a dynamic environment (e.g. use proportional weight defined by a weight-cycle-fraction, similar to the MCF described for development). Also to be considered is how to extend the moulting-weight concept to nauplii, whose weight ranges overlap among stages. Ultimately, use of moulting lengths may be a better approach, because length increases with stage even when weight is variable (van den Bosch and Gabriel, 1994). Errors in modelled development times are also caused in weight-based models by uncertainty in the growth model parameterizations. A seemingly trivial change in the allometric exponent is amplified and results in a 20% difference in modelled CV stage duration (e.g. 1 week at 48C/high food). Commonly used ranges of assimilation efficiency will alter modelled generation times by .2.5 weeks at 88C and high food, and longer for lower temperatures and food. Errors are compounded by

characterizations of environmental dependencies. Common parameterizations for temperature effects underestimate the true non-linearity, introducing errors of 30% into modelled stage durations and predicting generation times that are premature by 1 month at 48C. Formulations used for food effects greatly overestimate food limitation, such that, at 50 mg C l21, some modelled development times are 240% too long, and others incorrectly predict arrested development. Because weight-based models are highly sensitive to choices of parameter values, and modelled development is a by-product of modelled growth, it is difficult to calibrate these models for mean development times, let alone for development time variability. Part of the issue relates to the fact that growth and development are decoupled in a number of ways, including their different Q10 temperature dependencies (Campbell et al., 2001; Figure 2) and that development can proceed even when animals are starved (Crain and Miller, 2001). Hence, although they are useful tools for exploring the mechanistic processes controlling zooplankton growth, until we have a better quantitative understanding of growth processes and the relationship between growth and development, weight-based models are not ideal for simulating population-stage structure.

Limitations of stage-based population models using ODEs

Figure 11. Modelled development time (DTCVI) and standard deviation (sDTCVI) for new IBM compared with data for new stage-based IBM (in the section An alternative approach: stage-based models of individuals), with Dt ¼ 1 d. Data are for all five conditions shown in Table 1: 48C, 88C, and 128C at high food (H), and medium (M) and low (L) food at 88C.

Our analysis demonstrated that both the modelled mean development time and associated variability depend entirely on the structure of stage-based population models using ODEs (i.e. number of variables and parameter values). Accuracy requires the modelled stages to be subdivided into age classes and specification of the constant maturation rate (ri) as the reciprocal of age-class duration. Use of other published forms for ri results in errors in mean generation times of 2 –3 weeks at 88C and high food, and larger for lower temperatures and food. Although the models can be designed to match the data for mean development times, it is not possible to do the same for modelled variability. First, the mathematics yields smaller modelled C.V.DTi for later vs. earlier stages, which contrasts observations where C.V.DTCVI is typically larger than most other stages. Second, modelled variability depends on the fixed number of subdivisions (age classes and stages), which means that coefficients of variability cannot change with environmental conditions.

411

Modelling copepod development It may be possible to resolve these issues by using a large number of subdivisions such that modelled variability would be underestimated, and subsequently increasing the modelled variability to realistic levels by introducing a new “biodiffusivity” correlated with environmental factors. Another alternative would be to design an ODE model like those presented above in the note on hybrid model structures, where ageing is decoupled from moulting and moulting is distributed [i.e. ageing rate ¼ 1/ti, moulting rate ¼ Equation (11)]. However, both of these approaches would entail large increases in the number of modelled variables, which would increase the computational cost, especially when these models are coupled with physical transport. More practical alternatives may be the use of non-traditional ODEs, such as delay differential equations (McCauley et al., 1996) or framing ODE models in terms of mean-age (Hu et al., in press; see below).

Limitations of stage-based models using difference equations There are subtle issues related to the use of difference equations that can result in inaccurate modelled mean development times. First, the fraction that moults at a particular time should not be based on age-within-stage when the environmental conditions are dynamic; the MCF formulation should be employed. Second, even in a constant environment, the choice of model Dt affects predicted mean development times because of age discretization error. For a 13-stage model of C. finmarchicus, the discrete age error associated with Dt ¼ 1 d overestimates mean generation time by .1 week at 88C and high food, and this would be longer for lower temperatures and food. Reasonable estimates of mean times are only obtained when Dt , 0.1 d. Although a larger timestep may be used by keeping track of the fraction of Dt that is spent in the new stage after moult [Equation (10)], it is unclear how such a strategy should be implemented in applications accounting for individual variability. Larger time-steps may also be used in aggregated models for which the relative importance of age discretization error is reduced (Hu et al., in press), although aggregated stages are not appropriate for applications of stage-dependent phenologies or vertical distributions (Titelman and Fiksen, 2004). The requirement of a small time-step results in a large number of modelled variables introducing the same problems with computational cost as for the ODEs. It also results in a modelled variability far smaller than what is observed. Modelled variability can be improved by allowing distributed moulting (Figure 5b), although modellers must be careful to use the correct formula to avoid introducing errors into the mean development times, particularly when sSDi is high (e.g. generation times underestimated by 3 weeks at 88C/high food, and larger for lower temperature and food, Figure 9, Column c). However, use of probabilistic moulting functions can compound some of the other problems. For example, the number of modelled variables increases significantly because the model must track sufficient age classes to account for late developers. Use of total stage times that are 3 standard deviations greater than the mean and Dt ¼ 0.1 d requires 683 variables at 88C and high food. This requirement would be even greater for lower temperatures and food, and may still be insufficient because some individuals do not moult before reaching the last age class. Furthermore, there is a large increase in computational costs caused by both additional variables and the need to calculate cumulative probability functions for each age class of every stage at every time-step [Equation (12)].

These computational issues may be mitigated by a promising new “mean-age” approach (Hu et al., in press), which significantly reduces the number of modelled variables. The model of Hu et al. (in press) is an example of a stage-based population model using difference equations that tracks abundance in a single age class for each stage [Equation (7), mi ¼ 1 for all i, no j subscript is required]. Instead of describing the proportion that moults as a constant, as is traditionally done (i.e. fi ¼ Dt/SDi), they use a probabilistic fi that depends on an additional model variable: the mean age-within-stage. At each time-step, this mean age changes, owing to ageing of the animals that did not moult during the time-step and the addition of new recruits to the stage, whose age is zero. Thus, moulting is conceptually similar to the distributed moulting occurring over multiple age classes (Figure 5b), but mean age is continuous, not discrete. Accuracy of modelled mean times depends on how well the mean age represents the true distribution of age within stage. Hu et al. (in press) demonstrated that the mean-age approach resulted in accurate mean development times for Pseudocalanus spp., even with four aggregated stages (i.e. eight variables) and Dt ¼ 1 d. However, more study is needed to assess the ability of this approach to represent development-time variability accurately, as well as how cumulative probabilities should be altered for dynamic environments.

Advantages of our alternative: stage-based model of individuals Our model avoids the problems related to simulating development that arise with other traditional approaches. Modelled times do not depend on uncertain parameters or relationships, because the only parameters needing specification are the laboratory measurements of stage-duration means and variances. Environmental dependencies are consistent with observations de facto, making the model highly accurate for a range of constant laboratory conditions (Figure 11). Because it is formulated with MCF, it is appropriate for applications with dynamic environmental conditions; the only modification would be use of empirical relationships to describe the variation in stage duration distributions in a dynamic environment (e.g. the section Experimental data; Soussi et al., 1997). What is more, there are no time-step constraints for accurate mean times caused by age discretization (results shown in Figure 11 are with Dt ¼ 1 d), because it is possible to track the residual time after moult for each individual [Equation (10)]. Our approach is individual-based, allowing for simulation of individual history and statistical comparison of results with observed population distributions. It can be adaptable to any copepod species when similar laboratory data are available. The model is not limited to the use of the normal distribution, because Fitness can be used to assign individual stage durations from other cumulative probabilities distributions (e.g. gamma). Also, the model structure is ideal for coupling with weight-based models of individuals in a manner that keeps developmental timing accurate and independent of growth to allow for investigation of stage-dependent effects on growth (e.g. Kimmerer et al., 2007). IBMs can be coupled with physical transport in a Lagrangian sense (i.e. particle tracking), which is an intuitive and efficient approach that has been used successfully in a number of studies of copepod distributions (e.g. Hannah et al., 1998; Miller et al., 1998; Batchelder et al., 2002; Johnson et al.,

412 2006), and computations can be kept tractable by considering each particle as representative of more than one individual (Carlotti and Wolf, 1998; Miller et al., 1998). In addition, the designation of an individual’s stage-dependent Fitness provides an intuitive basis for characterizing other biological quantities (e.g. poorer Fitness may have higher mortality rates, etc.), as well as for investigations of the effects of an individual’s history or environmental factors on development. Therefore, the model presented here is an ideal tool that may be expanded to study prey consumption, dormancy, mortality rate estimation, naupliar development time, and how distributions of copepod abundance and production may be sensitive to climate change.

Acknowledgements We thank Gordon Fenton and Haiying Zhou for their guidance with the probabilistic formulae, and Morven Gentleman for his advice about parameter estimation and editorial feedback. We are grateful to Cabell Davis, Qiao Hu, Colleen Petrick, and two anonymous reviewers for their thoughtful comments that greatly improved the generality of the manuscript. We also acknowledge Holly Goulding, Jamie Pringle, Romero Advincula, and Chrissy Galloway for their input on early versions of the work. A National Science Foundation Grant to WCG (Biological Oceanography - 0222309) supported this research. This is GLOBEC contribution 586.

References Batchelder, H. P., Edwards, C. A., and Powell, T. M. 2002. 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. Progress in Oceanography, 53: 307 – 333. Batchelder, H. P., and Miller, C. B. 1989. Life history and population dynamics of Metridia pacifica: results from simulation modelling. Ecological Modelling, 48: 113 – 136. Batchelder, H. P., and Williams, R. 1995. Individual-based modelling of the population dynamics of Metridia lucens in the North Atlantic. ICES Journal of Marine Science, 52: 469 – 482. Campbell, R. G., Wagner, M. W., Teegarden, G. J., Boudreau, C. A., and Durbin, E. G. 2001. Growth and development rates of the copepod Calanus finmarchicus reared in the laboratory. Marine Ecology Progress Series, 221: 161– 183. Carlotti, F., Giske, J., and Werner, F. 2000. Modeling zooplankton dynamics. In ICES Zooplankton Methodology Manual, pp. 572– 667. Academic Press. Carlotti, F., Krause, M., and Radach, G. 1993. Growth and development of Calanus finmarchicus related to the influence of temperature: experimental results and conceptual model. Limnology and Oceanography, 38: 1125– 1134. Carlotti, F., and Nival, P. 1992. Model of copepod growth and development: moulting and mortality in relation to physiological processes during an individual moult cycle. Marine Ecology Progress Series, 84: 219 – 233. Carlotti, F., and Radach, G. 1996. Seasonal dynamics of phytoplankton and Calanus finmarchicus in the North Sea as revealed by a coupled one-dimensional model. Limnology and Oceanography, 41: 522– 539. Carlotti, F., and Sciandra, A. 1989. Population dynamics model of Euterpina acutifrons (Copepoda: Harpacticoida) coupling individual growth and larval development. Marine Ecology Progress Series, 56: 225 – 242. Carlotti, F., and Wolf, K. U. 1998. A Lagrangian ensemble model of Calanus finmarchicus coupled with a 1D ecosystem model. Fisheries Oceanography, 7: 191 – 204.

W. C. Gentleman et al. Caswell, H. 1989. Matrix Population Models: Construction, Analysis and Interpretation. Sinauer Associates, Sunderland, MA. Cook, K. B., Bunker, A., Hay, S., Hirst, A. G., and Speirs, D. C. 2007. Naupliar development times and survival of the copepods Calanus helgolandicus and Calanus finmarchicus in relation to food and temperature. Journal of Plankton Research, 29: 757– 767 Corkett, C. J., McLaren, I. A., and Sevigny, J. M. 1986. The rearing of the marine calanoid copepods Calanus finmarchicus (Gunnerus), C. Glacialis Jaschov, and C. Hyperboreous Krøyer with comment on the equiproportional rule. In Proceedings of the Second International Conference on Copepoda, 13 – 17 August 1984, pp. 539– 546 Ed. by G. Schriever, H. K. Schminke and C-T. Shih. National Museums of Canada, Ottawa. Crain, J. A., and Miller, C. B. 2001. Effects of starvation on intermolt development in a comparison between theoretical models and field studies. Deep Sea Research II, 48: 551 – 566. Davis, C. S. 1983. Laboratory rearing of marine calanoid copepods. Journal of Experimental Marine Biology and Ecology, 71: 119 – 133. Davis, C. S. 1984. Interaction of a copepod population with the mean circulation of Georges Bank. Journal of Marine Research, 42: 573– 590. Fiksen, Ø. 2000. The adaptive timing of diapause — a search for evolutionarily robust strategies in Calanus finmarchicus. ICES Journal of Marine Science, 57: 1825 – 1833. Gentleman, W., Leising, A., Frost, A., Strom, S., and Murray, J. 2003. Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics. Deep Sea Research Part II, 50: 2847 – 2875. Hannah, C. G., Naimie, C. E., Loder, J. W., and Werner, F. E. 1998. Upper-ocean transport mechanisms from the Gulf of Maine to Georges Bank, with implications for Calanus supply. Continental Shelf Research, 17: 1887– 1911 Heath, M., Robertson, W., Mardaljevic, J., and Gurney, W. S. G. 1997. Modelling the population dynamics of Calanus on the Fair Isle Current off northern Scotland. Journal of Sea Research, 38: 381– 412. Hind, A., Gurney, W. S., Heath, M., and Bryant, A. D. 2000. Overwintering strategies in Calanus finmarchicus. Marine Ecology Progress Series, 193: 95 – 107. Hirst, A. G., Peterson, W. T., and Rothery, P. 2005. Errors in juvenile copepod growth rate estimates are widespread: problems with the Moult Rate method. Marine Ecology Progress Series, 296: 263– 279. Hoaglin, D. C., Mosteller, F. M., and Tukey, J. W. 2000. Understanding Robust and Exploratory Data Analysis. John Wiley, New York. Hu, Q., Davis, C. S., and Petrik, C. M. A simplified age-stage model for copepod population dynamics. Marine Ecology Progress Series, in press. Hu, Q., Petrik, C. M., and Davis, C. S. 2007. Normal versus gamma: stochastic models of copepod molting rate. Journal of Plankton Research, 29: 985– 997. Johnson, C., Pringle, J., and Chen, C. 2006. Transport and retention of dormant copepods in the Gulf of Maine. Deep Sea Research II, 53: 2520– 2536. Kimmerer, W. J., Hirst, A. G., Hopcroft, R. R., and McKinnon, A. D. 2007. Estimating juvenile copepod growth rates: corrections, inter-comparisons and recommendations. Marine Ecology Progress Series, 336: 187–202. Klein Breteler, W. C. M., Schogt, N., and van der Meer, J. 1994. The duration of copepod life stages estimated from stage-frequency data. Journal of Plankton Research, 16: 1039 – 1057. Landry, M. R. 1975. The relationship between temperature and the development of life stages in the marine copepod Acartia clausi Giesbr. Limnology and Oceanography, 20: 854– 857. Leandro, S. M., Queiroga, H., Rodrigues-Grana, L., and Tiselius, P. 2006. Temperature-dependent development and somatic growth

Modelling copepod development in two allopatric populations of Acartia clausi (Copepoda: Calanoida). Marine Ecology Progress Series, 322: 189– 197. Lewis, C. V. W., Davis, C. S., and Gawarkiewicz, G. 1994. Wind forced biological-physical interactions on an isolated offshore bank. Deep Sea Research II, 41: 51 – 73. Lynch, D. R., Gentleman, W. C., McGillicuddy, D. J., and Davis, C. S. 1998. Biological/physical simulations of Calanus finmarchicus population dynamics in the Gulf of Maine. Marine Ecology Progress Series, 169: 189 – 210. Mathworks. 2005. MATLAB (Data Analysis Software System), Version 7.0. Natick, MA. McCauley, E., Nisbet, R. M., De Roos, A. M., Murdoch, W. W., and Gurnery, W. S. C. 1996. Structured population models of herbivorous zooplankton. Ecological Monographs, 66: 479 – 501. McLaren, I. A. 1986. Is “structural” growth of Calanus potentially exponential? Limnology and Oceanography, 31: 1342– 1346. McLaren, I. A. 1997. Modeling biases in estimating production from copepod cohorts. Limnology and Oceanography, 42: 584– 589. Miller, C. B., Johnson, J. K., and Heinle, D. R. 1977. Growth rules in the marine copepod genus Acartia. Limnology and Oceanography, 22: 326–335. Miller, C. B., Lynch, D. R., Carlotti, F., Gentleman, W., and Lewis, C. V. W. 1998. Coupling of an individual-based population dynamic model of Calanus finmarchicus to a circulation model for the Georges Bank region. Fisheries Oceanography, 7: 219 – 234. Miller, C. B., and Tande, K. S. 1993. Stage duration estimation for Calanus populations, a modelling study. Marine Ecology Progress Series, 102: 15 – 34. Peterson, W. T. 1986. Development, growth and survivorship of the copepod Calanus marshallae in the laboratory. Marine Ecology Progress Series, 29: 61 – 72. Peterson, W. T., and Painting, S. J. 1990. Developmental rates of the copepods, Calanus austalis and Calanoides carinatus in the laboratory, with discussion of methods used for calculation of development time. Journal of Plankton Research, 12: 283 – 293. Plaga´nyi, E. E., Hutchings, L., Field, J. G., and Verheye, H. M. 1999. A model of copepod population dynamics in the southern Benguella upwelling region. Journal of Plankton Research, 21: 1691– 1724. Sciandra, A. 1986. Study and modelling of the development of the copepod Euterpina acutifrons (Copepoda: Harpacticoida). Journal of Plankton Research, 8: 1149– 1162.

413 Soussi, S., and Ban, S. 2001. The consequences of individual variability in moulting probability and the aggregation of stages for modelling copepod population dynamics. Journal of Plankton Research, 23: 1279– 1296. Soussi, S., Carlotti, F., and Nival, P. 1997. Food- and temperaturedependent function of moulting rate in copepods: an example of parameterization for population dynamics models. Journal of Plankton Research, 19: 1331– 1346. Steele, J. H., and Frost, B. W. 1977. The structure of plankton communities. Philosophical Transactions of the Royal Society of London, Series B, Biological Sciences, 280: 485 –534. Steele, J. H., and Mullin, M. M. 1977. Zooplankton dynamics. In The Sea, vol. 6, pp. 857–890. Ed. by E. D. Goldberg. Interscience, New York. Titelman, J., and Fiksen, Ø. 2004. Ontogenetic vertical distribution patterns in small copepods: field observations and model predictions. Marine Ecology Progress Series, 284: 49 – 63. Tsuda, A. 1994. Starvation tolerance of a planktonic marine copepod Pseudocalanus newmani Frost. Journal of Experimental Marine Biology and Ecology, 181: 81 – 89. van den Bosch, F., and Gabriel, W. 1994. A model of growth and development in copepods. Limnology and Oceanography, 39: 1528– 1542. Vidal, J. 1980. Physioecology of zooplankton II. Effects of phytoplankton concentration, temperature, and body size on the development and molting rates of Calanus pacificus and Pseudocalanus sp. Marine Biology, 56: 185– 197. Wroblewski, J. S. 1980. A simulation of the distribution of Acartia clausi during Oregon Upwelling, August 1973. Journal of Plankton Research, 2: 43– 68. Wroblewski, J. S. 1982. Interaction of currents and vertical migration in maintaining Calanus marshallae in the Oregon upwelling zone — a simulation. Deep Sea Research, 29: 665– 686. Zakardjian, B. A., Sheng, J., Runge, J. A., McClaren, I., Plourde, S., Thompson, K. R., and Gratton, Y. 2003. Effects of temperature and circulation on the population dynamics of Calanus finmarchicus in the Gulf of St Lawrence and Scotian Shelf: study with a coupled, three-dimensional hydrodynamic, stage-based life history model. Journal of Geophysical Research, 108: 8016, doi:10.1029/2002JC001410. doi:10.1093/icesjms/fsn047