Copepod life strategy and population viability in response to prey timing and temperature: Testing a new model across latitude, time, and the size spectrum

A new model (“Coltrane”: Copepod Life-history Traits and Adaptation to Novel Environments) describes environmental controls on copepod populations via (1) phenology and life history and (2) temperature and energy budgets in a uniﬁed framework. The model tracks a cohort of copepods spawned on a given date using a set of coupled equations for structural and reserve biomass, developmental stage, and survivorship, similar to many other individual-based models. It then analyzes a family of cases varying spawning date over the year to produce population-level results, and families of cases varying one or more traits to produce community-level results. In an idealized global-scale testbed, the model correctly predicts life strategies in large Calanus spp. ranging from multiple generations per year to multiple years per generation. In a Bering Sea testbed, the model replicates the dramatic variability in the abundance of Calanus glacialis/marshallae observed between warm and cold years of the 2000s, and indicates that prey phenology linked to sea ice is a more important driver than temperature per se . In a Disko Bay, West Greenland testbed, the model predicts the viability of a spectrum of large-copepod strategies from income breeders with a adult size ∼ 100 µ gC reproducing once per year through capital breeders with an adult size > 1000 µ gC with a multiple-year life cycle. This spectrum corresponds closely to the observed life histories and physiology of local populations of Calanus ﬁnmarchicus , C. glacialis , and Calanus hyperboreus . Together, these complementary initial experiments demonstrate that many patterns in copepod community composition and productivity can be predicted from only a few key constraints on the individual energy budget: the total energy available in a given environment per year; the energy and time required to build an adult body; the metabolic and predation penalties for taking too long to reproduce; and the size and temperature dependence of the vital rates involved.

A new model ("Coltrane": Copepod Life-history Traits and Adaptation to Novel Environments) describes environmental controls on copepod populations via (1) phenology and life history and (2) temperature and energy budgets in a unified framework. The model tracks a cohort of copepods spawned on a given date using a set of coupled equations for structural and reserve biomass, developmental stage, and survivorship, similar to many other individual-based models. It then analyzes a family of cases varying spawning date over the year to produce population-level results, and families of cases varying one or more traits to produce community-level results. In an idealized globalscale testbed, the model correctly predicts life strategies in large Calanus spp. ranging from multiple generations per year to multiple years per generation. In a Bering Sea testbed, the model replicates the dramatic variability in the abundance of Calanus glacialis/marshallae observed between warm and cold years of the 2000s, and indicates that prey phenology linked to sea ice is a more important driver than temperature per se. In a Disko Bay, West Greenland testbed, the model predicts the viability of a spectrum of large-copepod strategies from income breeders with a adult size ∼100 µgC reproducing once per year through capital breeders with an adult size >1000 µgC with a multiple-year life cycle. This spectrum corresponds closely to the observed life histories and physiology of local populations of Calanus finmarchicus, C. glacialis, and Calanus hyperboreus. Together, these complementary initial experiments demonstrate that many patterns in copepod community composition and productivity can be predicted from only a few key constraints on the individual energy budget: the total energy available in a given environment per year; the energy and time required to build an adult body; the metabolic and predation penalties for taking too long to reproduce; and the size and temperature dependence of the vital rates involved.

INTRODUCTION
Calanoid copepods occupy a crucial position in marine food webs, the dominant mesozooplankton in many temperate and polar systems, important to packaging of microbial production in a form accessible to higher predators. They also represent the point at which biogeochemical processes, and numerical approaches like NPZ (nutrient-phytoplankton-zooplankton) models, start to be significantly modulated by life-history and behavioral constraints. The population-and community-level response of copepods to environmental change (temperature, prey availability, seasonality) thus forms a crucial filter lying between the biogeochemical impacts of climate change on primary production patterns and the food-web impacts that follow.
Across many scales in many systems, the response of fish, seabirds, and marine mammals to climate change has been observed, or hypothesized, to follow copepod community composition more closely than it follows total copepod or total zooplankton production. Examples include interannual variation in pollock recruitment in the Eastern Bering Sea Eisner et al., 2014), interdecadal fluctuations in salmon marine survival across the Northeast Pacific (Mantua et al., 1997;Hooff and Peterson, 2006;Burke et al., 2013), and long-term trends in forage fish and seabird abundance in the North Sea (Beaugrand and Kirby, 2010;MacDonald et al., 2015). These cases can be all be schematized as following the "junk food" hypothesis (Österblom et al., 2008) in which the crucial axis of variation is not between high and low total prey productivity, but rather between high and low relative abundance of large, lipid-rich prey taxa.
Calanoid copepods range in adult body size by more than two orders of magnitude, from <10 to >1000 µg C. Lipid content is likewise quite variable (Kattner and Hagen, 2009), even among congeneric species in a single environment (Swalethorp et al., 2011). Many but not all species enter a seasonal period of diapause in deep water, in which they do not feed and basal metabolism is reduced to ∼1/4 of what it is during active periods . Reproductive strategies include both income breeding (egg production fueled by ingestion of fresh prey during phytoplankton blooms) and capital breeding (egg production fueled by stored lipids in winter), as well as hybrids between the two strategies (Hirche and Kattner, 1993;Daase et al., 2013). Generation lengths vary from several weeks to several years.
These life-history traits (generation length, diapause, reproductive strategy, and annual routine more generally) constitute the mechanistic link between environment and the quality of the copepod community as prey (i.e., body size and composition). Lipid storage, coupled to diapause in deep water, is a strategy for surviving the winter in environments where winter foraging is not cost-effective energetically; and just as important, it provides energetic free scope for optimizing reproductive timing relative to prey availability (Falk-Petersen et al., 2009;Varpe et al., 2009). Lipid storage is tied to climate via temperature (which determines the rate at which an animal burns through its reserves during winter and rates of ingestion, growth, and development year-round) and phenology (i.e., timing of the copepods phytoplankton and protist prey: Mackas et al., 2012). This logic provides a route by which the energetics of fish, seabird, and mammal foraging are tied to temperature and phytoplankton phenology via the tradeoffs governing copepod life history.
There is likely a gap, then, between the focus of conventional oceanographic plankton models-total productivity by functional group-and the copepod traits of greatest importance to predators. A number of dynamical-modeling studies have attempted to fill this gap by modeling the copepods species by species in relation to climate forcing, often in an individualbased-model (IBM) framework (Miller et al., 2002;Ji et al., 2012;Maar et al., 2013;Wilson et al., 2016). There are two key limitations to the species-by-species approach, however. First, it is difficult to see how it can scale or generalize to the community level, given that our empirical information on the physiology and life history of the copepods is a patchwork, and realistically will always remain so. Second, it does not address the question of adaptation, either on the individual or species level. As individuals make use of their phenotypic plasticity in behavior, physiology, and life cycle, and as natural selection acts on existing species and subpopulations, it is likely that shifts in the biogeography of copepod traits such as size, lipid content, and life history pattern will not move in lockstep with the biogeography of existing species (Barton et al., 2013). Indeed, subpopulations of individual copepod species display so much life-history and physiological diversity (Heath et al., 2004;Daase et al., 2013) that it is not clear what the basic units of a general species-based model would even be. Observations of hybridization among species (Parent et al., 2015) only underscore this problem.
This paper presents a proof-of-concept for a trait-based, as opposed to species-based, copepod IBM, intended for eventual use in problems linking planktivores to climate and environment on global or regional scales. Record et al. (2013) presented a copepod community IBM in which explicit competition via a genetic algorithm was used to pick community assemblages out of a trait-based metacommunity along a latitudinal gradient. That study was concerned mainly with the emergent behavior of a very complex model system (predation-structured competition along with the interacting effects of six variable traits). In contrast, we have included as few explicitly variable traits as possible, guided by a strategic set of heuristic and quantitative comparisons with data (Figure 1). The balance point we have sought in this phase of work is the lightest-weight representation of diversity and plasticity that allows the model to (1) generate a realistic landscape of competitors in a single environment, (2) correctly predict fitness fluctuations in one population as a function of habitat, and (3) give sensible results over a wide biogeographic range.
The first of these criteria, captured by a Disko Bay, West Greenland model experiment (Figure 1, Section 3.4) is central to the goal of eventually allowing climate-to-copepod model studies to replace hand-picked sets of fixed types with a trait continuum. The second and third criteria (captured by a Bering Sea hindcast experiment and an heuristic, idealized biogeographic experiment: Figure 1 distinguish the effects of temperature and prey seasonality. As we will show, these initial experiments suggest a general hypothesis: that the viability of the calanoid community, at least near its highlatitude limit, is much more sensitive to prey abundance and phenology than to temperature.

General Approach
The model introduced here is "Coltrane" (Copepod Life-history Traits and Adaptation to New Environments) version 1.0. Matlab source code is available at https://github.com/neilbanas/coltrane. An overview of the model structure is shown in Figure 2. Like many individual-based models, Coltrane represents the time-evolution of one cohort of a clonal population, all bearing the same traits and spawned on the same date t 0 , with a set of ODEs. The state variables describing a cohort are relative developmental stage D, where D = 0 represents a newly spawned egg and D = 1 an adult; survivorship N, the fraction of initially spawned individuals that remain after some amount of cumulative predation mortality; structural biomass per individual S, and "potential" or "free scope" ϕ, which represents all net energy gain not committed to structure, or equivalently, the combination of internal energy reserves and eggs already produced. Combining reserves and eggs into one pool in this way lets us cleanly separate results that depend only on the fundamental energy budget (gain from ingestion, loss to metabolism, and energy required to build somatic structure) from results that depend on particular assumptions about egg production (costs, cues, and strategies). An alternate form of the model explicitly divides ϕ into internal reserves R and egg production rate E: the simpler model without this distinction will be called the "potential" or ϕ model and the fuller version the "egg/reserve" or ER model.
The ϕ and ER models take different approaches to generating population-level results from this cohort model, as explained in detail below (Section 2.4). In both cases, the logic changes from the simple forward time-integration at the cohort level: one runs the cohort model for all possible spawning dates t 0 , retroactively determines which spawning dates would prove optimal or sustainable, and considers the cohort time series from those t 0 values, appropriately weighted, to constitute the model solution (Section 2.4). The biological logic here is similar to the backwards-in-time dynamical optimization method frequently used in studies of optimal annual routines (Houston et al., 1993;Varpe et al., 2007), although our solving method is quite different and less exact. This is a compromise with the eventual goal of coupling Coltrane to oceanographic models as a spatially explicit IBM.
Communities are generated in Coltrane 1.0 simply by running families of cases of the population-level model that vary one or more traits. Treating coexisting populations as uncoupled vastly simplifies the interpretation of the landscape of viable strategies in a given environment, or the fundamental niche of a particular trait combination, our primary modes of analysis. At the same time, it tightly restricts our choices regarding the formulation of predation mortality. In reality, coupling through shared predators can rival bottom-up effects as a determinant of community structure (Holt et al., 1994;Chesson, 2000;Record et al., 2013Record et al., , 2014, and we expect that many potential applications of this model would require that this be better represented. In the present study, we have taken the minimalist, incrementalist approach of imposing the simplest possible form of predation mortality-a linear function, with scalings that closely mirror the growth and development functions (Section 2.2)-and restricting FIGURE 2 | Overview of the logic and execution of the two versions of the Coltrane 1.0 model. In the ϕ version, a cohort is represented by four state variables (D, S, ϕ, N) integrated forward in time over the cohort's lifespan. These are used to calculate a time series of fitness F. Next, this calculation is repeated across a family of spawning dates t 0 , and optimal and viable t 0 cases determined. This constitutes a population-level description, which then can be repeated across a range of trait values (u 0 , relative development rate) to describe a size-variable metacommunity. In the ER version, at the cohort level, ϕ is replaced by the state variable R and a time series of egg production E. Across a family of t 0 cases, a transition-matrix method is used to determine a stable annual pattern of relative egg production n(t 0 ), which is taken as the population-level prediction. A metacommunity is formed by varying two traits, u 0 and the date at which egg production begins t egg .
the terms of analysis. In particular, we will describe model output in terms of trait correlations, optimality, and viability, but not in terms of absolute copepod biomass or abundance. Likewise, while some plankton models resolve the process of adaptation explicitly (Clark et al., 2013), we address it only in the indirect sense of mapping the viable and optimal regions of the strategy landscape. This approach is less mechanistic but also helpfully agnostic about whether adaptation in the copepods arises through individual plasticity, species composition shifts, or natural selection per se.
An environment in Coltrane 1.0 is defined by annual cycles of three variables, total concentration of phytoplankton/microzooplankton prey P, surface temperature T 0 , and deep temperature T d . At present, these annual cycles are assumed to be perfectly repeatable, so that a "viable" strategy can be defined as a set of traits that lead to annual egg production above the replacement rate, given P, T 0 , and T d as functions of yearday t. The level of predation mortality (Section 2.2.4) might also be viewed as an environmental characteristic.

Ontogenetic Development
Calanoid copepods have a determinate developmental sequence, comprising the embryonic period, six naupliiar stages (N1-6), five copepodid stages (C1-5), and adulthood (C6). Similar to Maps et al. (2012), conversions between relative developmental stage D and the actual 13-stage sequence have been done using relative stage durations for C. finmarchicus from Campbell et al. (2001), which appear to be appropriate for other Calanus spp. with the proviso that C5 duration is particularly variable and strategy-dependent. Development in the model follows where developmental rate u is and All variables and parameters are defined in Table 1. Activity level a is, in this version of the model, a two-state switch calculated at each time step, 1 during active feeding and 0 during diapause. The temperature-dependent factor q d describes a power-law response with a Q 10 of Q d , where temperature is assumed to be T 0 during active feeding and T d during diapause. We use the Q 10 functional form for convenience: the differences between this and the leading alternatives (Belěhrádek, Arrhenius : Forster et al., 2011;Record et al., 2012) appear to be small compared with interspecies differences in this study . Prey saturation σ is a simple Michaelis-Menten function with half-saturation K s . The parameter u 0 , the development rate corrected to 0 • C, was found by  to be the primary trait responsible for differences in adult body size among Calanus spp. and other calanoids >50 µg C adult size, although not at a broader scale of diversity. It represents the aspect of development-rate variation that we interpret to be a strategy choice as opposed to a physiological or thermodynamic constraint.

Energy Gain and Loss
The two energy stores S (structure) and ϕ (reserves/potential) follow where G is net energy gain (ingestion minus metabolic losses). When net gain is positive, it is allocated between structure and potential according to the factor f s , which commits net gain entirely to structure before a developmental point D s , entirely to potential during adulthood, and to a combination of them in between: When G ≤ 0, the deficit is taken entirely from reserves: f s = 0. Before the first feeding stage (D < D f ) we assume G = 0 for simplicity. After feeding begins, where ingestion I and metabolic loss M are given by and r a is an assimilation efficiency. Ingestion follows a Kleiber's Law-like dependence on structural body mass S, with θ = 0.7 (Kleiber, 1932;Saiz and Calbet, 2007). I 0 is specific ingestion rate at saturating prey concentration, T = 0 • C, and S = 1 µg C. This is modulated by the activity switch a and prey saturation σ as in Equation (2), and a power-law temperature response for growth which is parallel to that for development (q d ) but with a different Q 10 . Q 10 values have been found to vary among copepod species but  argue that common values derived from a fit across community-level data are more appropriate for comparing species near their thermal optima. We use Q g = 2.5 and Q d = 3.0, as an approximation to the best-fit complex allometric curves reported by Forster et al. (2011). Energy loss to metabolism M follows the same temperature and size scalings. The factor r m is the ratio of metabolism to ingestion when prey is saturating. Unlike development and ingestion, which are assumed zero during diapause, M during diapause is nonzero but reduced to a basal fraction r b ≈ 1/4 : Note that in this formalism, gross growth efficiency ǫ becomes when a = 1. We have set r m = 0.14 such that ǫ = 0 when P = 1/4 K s .

Starvation
Potential ϕ is allowed to run modestly negative, to represent consumption of body structure during starvation conditions. A cohort is terminated by starvation if where in this study r starv = 0.1. A convenient numerical implementation of this scheme is to integrate S implicitly so that it is guaranteed > 0, and to integrate ϕ explicitly so that it is allowed to change sign, with no change of dynamics at ϕ = 0.

Predation Mortality
Predation mortality is assumed to have the same dependence on temperature and body size as ingestion, metabolism, and net gain (Hirst and Kiørboe, 2002). Survivorship N is set to 1 initially and decreases according to (it is convenient to calculate the numerical solution using ln N rather than N as the state variable, since values become extremely small). The mortality rate m is Frontiers in Marine Science | www.frontiersin.org 5 November 2016 | Volume 3 | Article 225 Banas and Campbell, 2016

K s
Half-saturation for ingestion See Table 2 m 0 Specific predation mortality at T = 0 • C, S = 1 µgC See such that that predation pressure relative to energy gain is encapsulated in a single parameter m 0 . In practice m 0 is a tuning parameter but we can solve for the value that would lead to an approximate equilibrium between growth and mortality. The condition is equivalent, by Equations (6) and (16), to and with a = 1 this becomes Averaging f s over the maturation period 0 ≤ D ≤ 1 with D s = 0.35, and assuming σ ≈ 2/3 on average for an organism that has aligned its development with the productive season, gives m 0 ≈ 0.2 I 0 . This is the default level of predation in the model except where otherwise specified.

Activity Level and Diapause
Modulation of activity level a has been treated as simply as possible, using a "myopic" criterion that considers only the instantaneous energy budget, rather than an optimization over the annual routine or lifetime (Sainmont et al., 2015). Furthermore, we treat a as a binary switch-diapause or full foraging activity-although intermediate overwintering states have been sometimes observed, e.g., C. glacialis/marshallae on the Eastern Bering Sea shelf in November (Campbell, personal communication). In the present model, we set a = 0 if D > D dia (the stage at which diapause first becomes possible) and the environment is such that total population energy gain would be higher under diapause. We can derive an expression for the threshhold at which this occurs by maximizing population energy gain as a function of a. When d/da of GSN is positive, active foraging a = 1 is the optimal instantaneous strategy and when it is negative, a = 0 is optimal. The threshhold can be rearranged to give a critical prey-saturation level where C dia = 1 + ϕ/S. The first term in Equation (23) can be derived more simply by setting dG/da = 0, a criterion based on ingestion and metabolism alone. The second term adjusts this criterion by discouraging foraging at marginal prey concentrations when predation is high. A third, temperaturedependent term has been neglected. The second, mortalitydependent term tends to produce unrealistic, rapid oscillations in which the copepods briefly "top up" on prey and then hide in a brief "diapause" to burn them. It is unclear whether this model behavior is a mathematical artifact-a limitation of combining actual lipid reserves and potential egg production into a single state variable-or whether it suggests that under some conditions the optimal level of foraging is intermediate between full activity and none. Incorporating a more mechanistic treatment of optimal foraging (Visser and Fiksen, 2013) and allowing a to vary continuously would address this. In this study, we have eliminated the phenomenon by approximating C dia as where r max ϕ = 1.5.

Eggs and Potential Eggs
The evolution equations above Equations (1), (6), (7), (16) specify the development of one cohort in the ϕ model. If this model is elaborated with an explicit scheme for calculating total egg production over time E(t), then it is possible to define R(t), individual storage/reserve biomass, and interpret R as a state variable and ϕ as a derived quantity. The relationship between the two is Thus, ϕ tracks the reserves that an animal would have remaining if it had not previously started egg production. This is a useful metric for optimizing reproductive timing, as we will show (Section 2.4). Any explicit expression for E(t) allows Equation (25) to replace Equation (7). In one model experiment below (Section 3.4), we use the following scheme: E(t) is the sum of income egg production E inc and capital egg production E cap , which are 0 until maturity is reached (D = 1) and an additional timing threshhold has been passed (t > t egg ). Past those threshholds, they are calculated as where E max is a maximum egg production rate which we assume to be equal to food-saturated assimilation: Thus, the trait t egg determines whether egg production begins immediately upon maturation (if t egg is prior to the date on which D reaches 1) or after some additional delay. Instead of t egg , expressed in terms of calendar day, one could introduce the same timing freedom through a trait linked to light, an ontogenetic clock that continues past D = 1, or a more subtle physiological scheme. However, since we run a complete spectrum of trait values in each environmental case, it is not important to the results how the delay is formulated, provided we only compare model output, rather than actual trait values, across cases.

Population-Level Response
A population-level simulation (Figure 2) consists of integrating either the ϕ model (Equations 1, 6, 7, 16) or ER model (Equations 1, 6, 16, 25) for a full annual cycle of spawning dates t 0 , and then identifying optimal and viable values of t 0 in terms of the egg fitness F, future egg production per egg (Varpe et al., 2007). Calculating a time series of F in the ϕ model requires an estimate of individual egg biomass W e in order to convert ϕ(t) from carbon units into a number of eggs, and a similar issue arises in the ER population model. Thus, a digression on the determination of W e is required.

Egg and Adult Size
The problem of estimating W e can be replaced by the problem of estimating adult size W a using the empirical relationship for broadcast spawners determined by Kiørboe and Sabatini (1995): where r ea = 0.013, θ ea = 0.62 (In the ER model, W a ≡ S + R at D = 1, but in the ϕ model we approximate it as S alone for simplicity). Adult size itself is an important trait for the model to predict, but the controls on it are rather buried in the model formulation above.  describe a theory relating body size to the ratio of development rate to growth rate based on a review of laboratory data for copepods with adult body sizes 0.3-2000 µgC. In our notation, their model can be derived as follows: if we approximate Equations (6), (7) in terms of a single biomass variable as then integrating from spawning to maturation gives since u is the reciprocal of the total development time. Growth rate has been written in terms of I 0 and an effective growth efficiency over the development period ǫ ′ . If we assume that egg biomass W e = W| D=0 is much smaller than W a = W| D=1 , then combining Equation (31) with Equation (2) gives Properly speaking, both ǫ ′ and T in Equation (32) are functions of t 0 since they depend on the alignment of the development period with the annual cycle. Since we are trying to use Equations (29), (32) to optimize t 0 , we have a circular problem. Record et al. (2013) derive an expression similar to Equation (32) and apply it iteratively because of this circularity. Some applications of Coltrane might require the same level of accuracy, but in the present study we take the expedient approach of simply assuming that T is the annual mean of T 0 and that ǫ ′ ≈ 1/3: i.e., that after t 0 is optimized, some diapause/spawning strategy will emerge that aligns the maturation period moderately well with a period of high prey availability. This assumption eliminates the need to run the model before estimating W e via Equations (29), (32).

Optimal Timing in the ϕ Model
With a method for approximating W e in hand, we can define egg fitness F as a function of ϕ. If a cohort spawned on t 0 were to convert all of its accumulated free scope ϕ-all net energy gain beyond that required to build an adult body structure-into eggs on a single day t 1 , the eggs produced per starting egg would be This expression condenses one copepod generation into a mapping F similar to the "circle map" of Gurney et al. (1992). Once the ODE model has been run for a family of t 0 cases, this mapping can be used to quickly identify optimal life cycles of any length. The optimal one-generation-per-year strategy is the t 0 that maximizes F 1 = F(t 0 → t 0 + 365). The optimal one-generation-per-two-years strategy has t 0 that maximizes F 1/2 = F(t 0 → t 0 + 2·365). The optimal two-generation-peryear strategy has spawning dates t 0 , t 1 that maximize the product F 2 = F(t 0 → t 1 ) · F(t 1 → t 0 + 365); and so on. A viable strategy is a combination of spawning dates and model parameters that give F ≥ 1.

Optimal Timing in the ER Model
In reality, of course, copepods are not free to physically store indefinite amounts of reserves within their bodies and then instantaneously convert them into eggs when the timing is optimal. If a scheme for calculating egg production over time E(t) is added to the model as in Section 2.3, then the per-generation mapping represented by F takes a different form. First, for each cohort, we use the assumption that the environmental annual cycle repeats indefinitely to convert the time series of ENegg production discounted by survivorship-to a function of yearday, by adding the value on days 365+i, 2·365+i, ... to the value on day i (in practice we discretize the year into 5 d segments rather than yeardays per se). Next, we construct a matrix V whose rows are the year-long time series of EN/W e for each spawning date t 0 . V is thus a transition matrix with spawning date in generation k running down rows and spawning date in generation k+1 running across columns. Given a discrete annual cycle n k of eggs spawned in generation k, one can calculate the expected annual cycle of egg production in the next generation as n k + 1 = V · n k , where n is given as a column vector. The first eigenvector of V gives a seasonal pattern of egg production that is stable in shape, with the corresponding eigenvalue λ giving one plus the population growth rate per generation: A strict criterion for strategy viability would then be λ ≥ 1, although this criterion is unhelpfully sensitive to predation mortality. A more robust criterion (which we use in Section 3.4 below) is to consider a strategy viable if it yields lifetime egg production above the replacement rate: if E(t 0 ; t) and N(t 0 ; t) are the time series of egg production and survivorship for a cohort spawned on t 0 , and n(t 0 ) is a normalized annual cycle of egg production, Thus, in the ER version of the model, as in the ϕ version, we have an efficient method that describes the long-term viability of a trait combination under a stable annual cycle, along with the optimal spawning timing associated with those traits in that environment; and these methods only require us to explicitly simulate one generation.

Assembling Communities
Community-level predictions in Coltrane take the form of bounds on combinations of traits that lead to viable populations in a given environment (Figure 2). There are many copepod traits represented in the model that one might consider to be axes of diversity or degrees of freedom in life strategy: u 0 , I 0 , θ , D s , K s , W e /W a , and even m 0 to the extent that predation pressure is a function of behavior (Visser et al., 2008). Record et al. (2013) allowed five traits to vary among competitors in their copepod community model. We have taken a minimalist approach, where in the ϕ model we allow only one degree of freedom: variation in u 0 from 0.005 to 0.01 d −1 .  showed from a review of lab studies that u 0 variations appear to be the primary mode of variation in adult size among large calanoids (W a > 50 µgC) including Calanus and Neocalanus spp., with slower development leading to larger adult sizes. That study also suggests that variation in I 0 is responsible for copepod size diversity on a broader size or taxonomic scale (e.g., between Calanus and small cyclopoids like Oithona). However, variation in I 0 (energy gain from foraging) probably only makes sense as part of a tradeoff with predation risk or egg survivorship (Kiørboe and Sabatini, 1995) and we have left the formulation of that tradeoff for future work. We therefore expect Coltrane 1.0 to generate analogs for large and small Calanus spp. (∼100-1000 µgC adult size) but not analogs for Oithona spp. or even small calanoids like Pseudocalanus or Acartia. Choices regarding reproductive strategy require another degree of freedom. In the ϕ model, this does not require additional parameters, because the difference between, e.g., capital spawning in winter and income spawning in spring is simply a matter of the time t at which F is evaluated in postprocessing: each model run effectively includes all timing possibilities (Equation 33). In the ER model, however, diversity in reproductive timing must be made explicit. Under the simple scheme for egg production specified above (Section 2.3), this takes the form of running a family of cases varying t egg for each t 0 and u 0 .

Model Experiments
This study comprises three complementary experiments ( Table 2). The first of these is an idealized global testbed which addresses broad biogeographic patterns. The second is a testbed representing the Eastern Bering Sea shelf, which addresses time-variability in one population in one environment. The last is a testbed representing Disko Bay, West Greenland, which addresses trait relationships along the size spectrum in detail. The first two are evaluated entirely in terms of the ϕ model, while in the Disko Bay case we use the ER model to allow more specific comparisons with observations.
The global testbed consists of a family of idealized environments in which surface temperature T 0 is held constant, and prey availability is a Gaussian window of width δt ′ centered on yearday 365/2: We compare environmental cases in terms of T 0 and an effective season length which rescales the δt ′ cases in terms of the equivalent number of days of saturating prey per year. We assume that deep, overwintering temperature T d = 0.4 T 0 . The ratio 0.4 matches results of a regression between mean temperature at 0 and 1000 m in the Atlantic between 20 and 90 • N, or 0 and 500 m in the Northeast Pacific over the same latitudes (World Ocean Atlas 2013: http://www.nodc.noaa.gov/OC5/woa13/). Over the same data compilation, the mean seasonal range in temperature is approximately 5 • C at the surface (and approximately zero at 500-1000 m); an alternate formulation of the testbed that models T 0 as an annual sinusoid with this range gives results that are somewhat noisier but heuristically very similar to those shown in Section 3.2 below. The Bering Sea testbed considers interannual variation in temperature, ice cover, and the effect of ice cover on in-ice and pelagic phytoplankton production (Stabeno et al., 2012b;Sigler et al., 2014;. Variation between warm, low-ice years and cold, high-ice years has previously been linked to the relative abundance of large zooplankton including C. glacialis/marshallae , and we test Coltrane predictions against 8 years of C. glacialis/marshallae observations from the BASIS program. Seasonal cycles of T 0 , T d , and P are parameterized using empirical relationships between ice and phytoplankton from Sigler et al. (2014) and a 42-year physical The Disko Bay testbed represents one seasonal cycle of temperature and phytoplankton and microzooplankton prey, based on the 1996-1997 time series described by Madsen et al. (2001). We use this particular dataset not primarily as a guide to the current or future state of Disko Bay but rather as a specific circumstance in which the life-history patterns of three coexisting Calanus spp. (C. finmarchicus, C. glacialis, C. hyperboreus) were documented (Madsen et al., 2001). Details are given in Section 3.4 and the Appendix in Supplementary Material.

An Example Population
One case from the global experiment with u 0 = 0.007 d −1 , T 0 = 1 • C, and δt = 135 is shown in detail in Figure 3 to illustrate the analysis method described in Section 2.4.2. In this case, out of cohorts spawned over the full first year, only those spawned in spring reached adulthood without starving ( Figure 3B, bluegreen lines; non-viable cohorts not shown). The fitness function F (Equation 33) declines during winter diapause and rises during the following summer when prey are available. There is no equivalent peak during the third summer, indicating that by this time cumulative predation mortality is so high that there is no net advantage to continuing to forage before spawning.
The maximum value of F for most cohorts ( * , Figure 3C) comes at ∼1.5 year into the simulation, at the peak in prey availability following maturation. This point in the annual cycle, however, does not fall within the window of spawning dates at which maturation is possible (compare year 2 in Figure 3C with year 1 in Figure 3B), and thus is an example of "internal life history mismatch" (Varpe et al., 2007), the common situation in which the spawning timing that maximizes egg production by the parent is not optimal for the offspring. The long-term egg fitness corresponding to stable 1-year and 2-year cycles is marked for each cohort (Figure 3C, red, orange circles). Some but not all of the cohorts that reach maturity are able to achieve F > 1, egg production above the replacement rate, in these cyclical solutions (solid circles). The best 1-year and 2-year strategies achieve similar maximum fitness values (red vs. orange solid dots), although they require slightly different seasonal timing.
Note that although F can be described as the egg-fitness function, the lines in Figure 3C-time series of F for particular spawning dates t 0 -are not the same as the seasonal curve of egg fitness that results from a backwards-in-time dynamic optimization (e.g., Figure 6F in Varpe et al., 2007). Rather, each curve of F(t 0 ; t) in our approach gives a series of possible values for egg fitness at t 0 depending on what future strategy is taken. The forwards and backwards calculations converge (at least qualitatively) once the internal life-history mismatch is resolved and a stable long-term cycle is found (red and orange circles, Figure 3C). As expected (Varpe et al., 2007), these stable values of egg fitness peak, for each generation length, somewhat prior to the bloom maximum ( Figure 3C).

Global Behavior
In the global experiment, populations like that shown in Figure 3 were run for a spectrum of u 0 values, across combinations of T 0 and δt from −2 to 16 • C and 0 to 310 d (the latter corresponding to δt ′ from 0 to 150 d). Across these cases, at a given u 0 , the model predicts a log-linear relationship between adult size and temperature, which is not much perturbed by variation in prey availability (Figure 4). The slope of this relationship is equivalent to a Q 10 of 1.8-2.0, consistent with that predicted by Equation (32): Field observations of size in relation to temperature in C. finmarchicus and C. helgolandicus across the North Atlantic show a similar relationship (Q 10 = 1.65, Wilson et al., 2015, with prosome length converted to carbon weight based on Runge et al., 2006). Somewhat surprisingly, even wide variation in prey conditions (clusters of gray dots, Figure 4) has only minor effects on this slope. The intercept of the size-temperature relationship depends on u 0 (Figure 4), with u 0 = 0.005-0.01 d −1 corresponding to the Developmental stage D for cohorts that reach maturity (D = 1) without starving: the blue-yellow color scale corresponds to spawning dates t 0 over the viable period from pre-bloom to bloom maximum. (C) Fitness F over time for the cohorts shown in (B), i.e., the expected value (eggs egg −1 ) of converting all free scope ϕ to eggs on a given date. Curves of F begin when maturity is reached (D = 1) and egg production becomes possible. Open and solid circles mark the value of F on the 1-year (orange) and 2-year (red) anniversaries of the original spawning date. Solid circles mark cohorts that achieve a fitness above the replacement rate. range of adult size from C. finmarchicus to C. hyperboreus at the cold end of the temperature spectrum (Disko Bay, ∼ 0 • C: Swalethorp et al., 2011). It is not always fair, however, to associate a particular u 0 value with a particular species over the full range of temperatures included. As  discuss further, the temperature response of an individual species is often dome-shaped, a window of habitat tolerance (Møller et al., 2012;Alcaraz et al., 2014), whereas Coltrane 1.0 uses the monotonic, power-law response observable at the community level (Forster et al., 2011). C. finmarchicus, for example, is fit well by u 0 = 0.007 d −1 at higher temperatures (4-12 • C), whereas near 0 • C in Disko Bay, it has been observed to be considerably smaller than extrapolation along the u 0 = 0.007 d −1 power law would predict. Past studies have also found C. finmarchicus growth and ingestion FIGURE 4 | Relationship between adult size W a and mean surface temperature T 0 in the "global" model experiment, for three values of relative development rate u 0 , in comparison with observations (Peterson, 1986;Swalethorp et al., 2011;Wilson et al., 2015;Campbell et al., in press) and laboratory results (Campbell et al., 2001;Rey-Rassat et al., 2002). Model results (gray dots) represent structural biomass S, as do observations marked with a ⋆; observations marked with a • represent total biomass R + S. Clusters of gray dots indicate families of model cases varying productive season length (horizontal axis in Figure 5).
to be suppressed at low temperatures, i.e., to show a very high Q 10 compared with the community-level value (Campbell et al., 2001;Møller et al., 2012).
With this caveat on the interpretation of u 0 , we can observe a sensible gradation in life strategy along the u 0 axis (Figure 5). From u 0 = 0.01 d −1 (C. finmarchicus-like at 0 • C) to u 0 = 0.005 d −1 (C. hyperboreus-like), the environmental window in which multi-year life cycles are viable (F 1/2 ≥ 1) expands dramatically. This window overlaps significantly with the window of viability for 1-year life cycles (F 1 ≥ 1; Figure 5, black vs. gray contours). In all u 0 cases, there is a non-monotonic pattern in maximum fitness as a function of either temperature or prey (Figure 5, color contours), as environments align and misalign with integer numbers of generations per year or years per generation.
The overall gradient from high to moderate F with increasing temperature (Figure 5) is largely an artifact of displaying F normalized to generation as opposed to per calendar year. In general, these results should not be taken as a quantitative prediction of annual production rates: the linear mortality closure that simplifies the analysis also omits the role of density dependence in stabilizing growth rates. Accordingly, in what follows, we consider only whether F in a given circumstance exceeds replacement rate, not whether it exceeds it modestly or dramatically.
The number of generations per year in the timing strategy that optimizes F for each (T 0 , δt) habitat combination is shown in Figure 6 for u 0 = 0.007 d −1 . This u 0 value corresponds in adult size to Arctic C. glacialis and temperate C. marshallae populations in the Pacific (Figure 4), species which coexist and are nearly indistinguishable in the Bering Sea. In the lowest-prey conditions, no timing strategy is found to be viable. As prey and temperature increase, the model predicts bands proceeding monotonically from multiple years per generation to multiple generations per year. Validating these model predictions requires parameterizing places (in terms of T 0 and δt) in addition to parameterizing their inhabitants, and thus the meaning of either success of failure is ambiguous. Still, we can observe the following. Ice Station Sheba in the high Pacific Arctic (Figure 1) falls in the non-viable regime (Figure 6), consistent with the conclusion of   Ashjian et al. (2003) that Calanus spp. are unable to complete their life cycle there. Disko Bay falls on the boundary of 1and 2-year generation lengths, consistent with observations of C. glacialis there (Madsen et al., 2001). At Newport, Oregon, near the southern end of the range of C. marshallae, the model predicts multiple generations per year, consistent with observations by Peterson (1979).

A High-Latitude Habitat Limit in Detail: The Eastern Bering Sea
These idealized experiments (Figures 5, 6) suggest that very short productive seasons place a hard limit on the viability of Calanus spp., regardless of size, temperature, generation length, or match/mismatch considerations (although these factors affect where exactly the limit falls). A decade of observations in the Eastern Bering Sea provide a unique opportunity to resolve this viability limit with greater precision. This analysis takes advantage of the natural variability on the Southeastern Bering Sea shelf described by the "oscillating control hypothesis" of Hunt et al. (2002Hunt et al. ( , 2011: in warm, low-ice years, the spring bloom in this region is late (∼ yearday 150; Sigler et al., 2014) and the abundance of large crustacean zooplankton including C. glacialis/marshallae is very low, while in colder years with greater ice cover, the pelagic spring bloom is earlier, ice algae are present in late winter, and large crustacean zooplankton are much more abundant. The task of replicating these observations serves to test the Coltrane parameterization, and situating them within a complete spectrum of temperature/ice cover cases also allows the model to provide some insight into mechanisms.
Mean surface temperature T 0 was used to index annual cycles of surface and bottom temperature on the Eastern Bering Sea middle shelf (Appendix in Supplementary Material; insets in Figure 7). Date of ice retreat t ice was likewise used to index phytoplankton availability over each calendar year (Appendix in Supplementary Material; insets in Figure 7). Coltrane was run for each (T 0 , t ice ) combination with u 0 = 0.007 d −1 , thus consistent with Figure 6 except for the more refined treatment of Frontiers in Marine Science | www.frontiersin.org 13 November 2016 | Volume 3 | Article 225 FIGURE 7 | Results of the "Bering" experiment. Color contours give the predicted egg fitness of a C. glacialis/marshallae analog under combinations of ice retreat timing (assumed to control spring bloom timing: Appendix in Supplementary Material) and temperature. Examples of the annual cycles of prey availability P and surface and bottom temperature T 0 , T d are given at left and bottom. Dots locate years 1971-2012 in this timing/temperature parameter space, for the northern (blue) and southern (red) middle shelf. Numbers beside red dots give the measured summer abundance of C. glacialis/marshallae on the southern shelf across years. White contours give the bounds of the central region within which Coltrane predicts one generation per year to be the optimal generation length.
environmental forcing, and an adjustment to K s to match results of Bering Sea feeding experiments (Campbell et al., in press). The maximum egg fitness F for a one-generation-per-year strategy is shown as a function of T 0 and t ice in the main panel of Figure 7. Coltrane predicts that one generation per year is the optimal life cycle length everywhere in this parameter space except for the cold/ice-free and warm/high-ice-cover extremes (white contours), combinations which do not occur anywhere in a model hindcast of middle-shelf conditions back to 1971 (Figure 7, red and blue dots). Late summer measurements of C. glacialis/marshallae abundance (individuals m −2 ), averaged over the middle/outer shelf south of 60 • N, are shown in Figure 7 for 2003-2010 (n = 364 over the 8 years; Eisner et al., 2014Eisner et al., , 2015. Both these observations and the predicted maximum F from Coltrane show a dramatic contrast between the warm years of 2003-05 (t ice = 0) and the cold years of 2007-2010 (t ice = 100-130), with the transitional year 2006 harder to interpret. Eisner et al. (2014) found that there was less contrast between cold year/warm year abundance patterns on the northern middle/outer shelf, consistent with the model prediction that all hindcast years on the northern shelf fall within the "viable" habitat range for C. glacialis/marshallae (Figure 7, blue dots).
The viability threshhold that the Southeastern Bering Sea appears to straddle is qualitatively similar to that in the more idealized global experiment (Figures 5, 6), primarily aligned with the phenological index (horizontal axis) rather than the temperature index (vertical axis). The threshhold in the Bering Sea experiment (t ice ≈ 90-100) falls somewhat beyond the dividing line imposed in the experiment setup between early, ice-retreat-associated blooms and late, open-water blooms (t ice = 75: see Appendix in Supplementary Material, Sigler et al., 2014). This gap (whose width depends on the mortality level m 0 : not shown) indicates that some period of ice algae availability is required by C. glacialis/marshallae in this system, in addition to a favorable pelagic bloom timing.

Coexisting Life Strategies in Detail: Disko Bay
The experiments above test the ability of Coltrane 1.0 to reproduce first-order patterns in latitude and time but do not provide sensitive tests of the model biology. A model case study in Disko Bay, where populations of three Calanus spp. coexist and have been described in detail (Madsen et al., 2001;Swalethorp et al., 2011), allows a closer examination of the relationships among traits within the family of viable life strategies predicted by Coltrane.
The model forcing (Figure 8) describes a single annual cycle, starting with the 1996 spring bloom. This represents a cold, highice state of the system, compared with more recent years in which the spring bloom is earlier (e.g., 2008, Figure 8, Swalethorp et al., 2011) and the deep layer is warmed by Atlantic water intrusions (Hansen et al., 2012). This particular year was chosen because measurements of prey availability and Calanus response by Madsen et al. (2001) were particularly complete and coordinated. A simple attempt to correct the prey field for quality and Calanus preference was made by keeping only the >11 µm size fraction of phytoplankton and adding total microzooplankton, in µg C. The measured phytoplankton C:chl ratio was used to convert the sum to an equivalent chlorophyll concentration, and this time series was then slightly idealized for clarity (Figure 8, Appendix in Supplementary Material).
Sensible results were only possible after tuning the predation mortality scale coefficient m 0 . It is likely that our simple mortality scheme introduces some form of bias, compared with the reality in this system of predation by successive waves of visual and non-visual predators, which will be considered in a separate study. Still, a sensitivity experiment using the ϕ model shows that varying m 0 has, as intended, a simple, uniform effect on fitness/population growth (Figure 9) that leaves other trait relationships along the size spectrum unaffected. The ϕ model predicts that copepods similar to C. finmarchicus in size have much greater fitness at a generation length of 1 year than at 2 years or more; that C. hyperboreus would be unable to complete its life cycle in 1 year, but is wellsuited to a 2-year cycle; and that C. glacialis falls in the size range where 1-and 2-year life cycles have comparable fitness value. These results are consistent with observations (Madsen et al., 2001) and more general surveys of life strategies in the three species (Falk-Petersen et al., 2009;Daase et al., 2013).  These results naturally raise the question of whether even lower values of u 0 -further reductions in development ratewould produce even larger copepods with even longer life cycles in this environment. C. hyperboreus has been reported to have a life cycle of up to 5 years in other systems (Falk-Petersen et al., 2009) and so the question is more than theoretical. In this version of Coltrane, the lower limit on development rate (and thus the upper limit on adult size) are set by the assumption that modulation of this rate is spread uniformly across the developmental period, rather than concentrated in late copepodid stages, as might be more realistic (Campbell et al., 2001). The largest viable adults in the Disko experiment are those that barely reach D = D s , the start of reserve accumulation, before a first diapause is required.
For greater specificity, we switched from the ϕ to the ER model version, running a spectrum of t egg cases (earliest possible eggproduction date: see Section 2.5) along with a spectrum of u 0 (development rate) cases. The predicted "community, " then is the set of all combinations of u 0 and t egg that lead to a viable level of lifetime egg production. An analog for each of the three Calanus spp. is constructed by averaging model results over the set of viable (u 0 , t egg ) cases that predict an adult size within 30% of the average measured adult size for that species. The ER model imposes additional constraints on the model organisms-e.g., they are no longer allowed an infinite egg production rate-and to compensate we reduced m 0 from 0.08 d −1 to 0.06 d −1 .
FIGURE 10 | Emergent relationship between generation length and adult size in the Disko Bay model experiment. Large colored dots indicate results for trait combinations that achieve a viable rate of egg production per generation (color coding matches that in Figure 12: blue, generation length of 1 year; light purple, 2 years; dark purple, 3 years) while small gray dots indicate trait combinations that reach maturity without starvation but have egg production rates below replacement level.
The relationship between generation length and adult size across all (u 0 , t egg ) combinations is shown in Figure 10. Results are consistent with the ϕ model (Figure 9) only a 1-year life cycle is viable for C. finmarchicus in this environment, only a 2-year or longer cycle is viable for C. hyperboreus, and C. glacialis again lies near the boundary where the two strategies are comparable. Note that the model allows for a continuum of intermediate cases in the C. finmarchicus-C. glacialis size range, consistent with the observation of hybridization between these species (Parent et al., 2015).
The ER model also predicts a time series of egg production associated with each trait combination, which we can compare with observations for each species. The model predicts that C. finmarchicus analogs spawn in close association with the spring bloom, that C. hyperboreus spawns well before the spring bloom, and that C. glacialis is intermediate (Figures 11,  12A). These patterns are all in accordance with Disko Bay observations (Madsen et al., 2001;Swalethorp et al., 2011), although the absolute range is muted: Madsen et al. (2001) report C. hyperboreus spawning as early as February. As one would expect from these timing patterns, the model predicts a significant trend between size and the capital fraction of total egg production E cap /(E inc + E cap ) ( Figure 12C). Again, the pattern is qualitatively correct but muted: Coltrane predicts 80% income breeding at the size of C. finmarchicus (a pure income breeder in reality) and 80% capital breeding at the size of C. hyperboreus (a pure capital breeder in reality). More notable than the error is how much of the income/capital spectrum can apparently be reproduced as a consequence of optimizing reproductive timing alone Varpe et al. (2009), without imposing the physiological difference between the two strategies as an independent trait (Ejsmond et al., 2015).
The model predicts ( Figure 12B) that the largest model organisms, with the longest generation lengths, enter their first diapause near the boundary between copepodite stages C4 and C5 (D ≈ 0.75), whereas smaller organisms enter first diapause well into stage C5. Madsen et al. (2001) found that both C. glacialis and C. hyperboreus diapause as C4, C5, and adults in Disko Bay, suggesting that the model is biased toward fast maturation. The discrepancy could also be related to intraspecific variation in the real populations or nonequiproportional development in the late stages, i.e., a variable conversion scale between actual developmental stage and D.
Finally, the ER version of Coltrane allows an estimate of the fraction of individual carbon in the form of storage lipids R/(R + S) ( Figure 12D). Averaging each model population from the first diapause-capable stage D = D dia through adulthood, weighted by survivorship N, yields an overall range that compares well with the species-mean wax ester fractions measured by Swalethorp et al. (2011): ∼30% for C. finmarchicus to ∼60% for C. hyperboreus. In the middle of the size spectrum, reserve fraction is highly variable across viable 2-year strategies, a warning that the success of this final model prediction may be partly fortuitous. Still, taken as a whole, this experiment has yielded a striking result: that a small set of energetic and timing contraints is able to correctly predict, a priori, that Disko Bay should be able to support a spectrum of calanoid copepods from income breeders with an adult size ∼100 µg C, a 1-year life cycle, and a wax ester fraction ∼30% to capital breeders with an adult size ∼1000 µg C, a two-or-more-year life cycle, and a wax ester fraction ∼60%.

Temperature and Timing
In the results above, whether prey availability is treated simply (Figures 5, 6) or with site-specific detail (Figure 7), it appears that the viability of the calanoid community near its highlatitude limit is more sensitive to prey abundance and phenology than to temperature. This result is heuristically similar to the conclusions of Ji et al. (2012) and Feng et al. (2016), although the exact physiological mechanisms differ. Alcaraz et al. (2014) suggested based on lab experiments that C. glacialis reaches an bioenergetic limit near 6 • C, and Holding et al. (2013) and others have hypothesized that thermal limits will produce ecosystemlevel tipping points in the warming Arctic. Our results, in contrast, suggest that thermal tipping points, even if present at the population level, do not generalize to the community level in copepods. Rather, the model predicts complete continuity between the life strategy of Arctic C. glacialis and temperate congeners like C. marshallae (Figure 6). It also suggests that even on the population level in the Bering Sea, warm/coldyear variation in prey availability is a sufficient explanation of variability in the abundance of C. glacialis/marshallae (Figure 7), without the invocation of a thermal threshhold.
Both the global and Bering experiments suggest, furthermore, that increasing water temperature per se is not necessarily a stressor on copepod communities, even high-latitude communities. In both cases, the low-prey viability threshhold actually relaxes (i.e., is tilted toward lower prey values) as temperature increases, indicating that in these testbeds, the positive effect of temperature on growth and maturation rate actually outweighs the effect of temperature on metabolic losses and overwinter survival (This result may be reliant on the model assumption that the stage of first diapause is highly plastic). In cases where deep, overwintering temperatures increase faster than surface temperatures (Hansen et al., 2012) this balance may not hold, and in the real ocean changes in temperature are highly confounded with changes in phytoplankton production and phenology. Still, it is notable that the model predicts that warming temperatures will have a non-monotonic effect on copepod populations (∂F/∂T 0 ≷ 0, Figures 5, 6) even when metabolic thermal threshholds sensu Alcaraz et al. (2014) and changes in prey availability are not considered. These results are a caution against overly simple climate-impacts projections based on temperature alone.

Uncertainties and Unresolved Processes
The biology in Coltrane could be refined in many ways, but two issues stand out as being both mechanistically uncertain Frontiers in Marine Science | www.frontiersin.org and sensitive controls on model behavior. These correspond to the two parameters that it was necessary to tune among model experiments ( Table 2) the obstacles to formulation of a fully portable scheme that could produce accurate results across the full range of environments considered here with a single parameterization.
The first of these is the perennial problem of the mortality closure. We modeled predation mortality as size-dependent according to the same power law used for ingestion and metabolism, a choice which is mathematically convenient and makes the effect of top-down controls, if not minor, then at least simple and easy to detect (Figure 9). This size scaling is consistent with the review by Hirst and Kiørboe (2002) but that study also shows that the variation in copepod mortality not explained by allometry spans orders of magnitude (cf. Ohman et al., 2004). Indeed, in some cases one might posit exactly the opposite pattern, in which mortality due to visual predators like larval fish increases with prey body size (Fiksen et al., 1998;Varpe et al., 2015). This latter pattern is one hypothesis for why in reality C. hyperboreus is confined to high latitudes, whereas the model predicts no southern (warm, high-prey) habitat limit to C. hyperboreus analogs based on bottom-up considerations (Figure 5). Merging Coltrane 1.0 with a light-and size-based predation scheme similar to Varpe et al. (2015) or Ohman and Romagnan (2015) would allow one to better test the balance of bottom-up and top-down controls on calanoid biogeography.
Second, our experience constructing the Bering Sea and Disko Bay cases suggests that the greatest uncertainty in the model bioenergetics is actually not the physiology itselfempirical reviews like Saiz and Calbet (2007), Maps et al. (2014), Kiørboe and Hirst (2014), and  have constrained the key rates moderately well-but rather the problem of translating a prey field into a rate of ingestion. Within each of our model testbeds, the prey time series P remains subject to uncertainty in relative grazing rates on ice algae, large and small pelagic phytoplankton, and microzooplankton, despite a wealth of local observations and a history of work on this problem in Calanus specifically (Olson et al., 2006;Campbell et al., 2009, in press). The precision of each testbed, and even moreso the ambition of generalizing across them, is also limited by uncertainty in the shape of the functional response (Frost, 1980;Gentleman et al., 2003), here represented by a halfsaturation coefficient, which has not been found to be consistent across site-specific studies (Campbell et al., in press;Møller et al., 2016) or well-constrained by general reviews . This ambiguity is perhaps not surprising when one considers that ingestion as a function of chlorophyll or prey carbon is not a simple biomechanical property, but in fact a plastic behavioral choice. Accordingly, it might well be responsive not only to mean or maximum prey concentration but also to the prey distribution over the water column, the tradeoff between energy gain and predation risk (Visser and Fiksen, 2013), prey composition and nutritional value, and the context of the annual routine. These issues are fundamental to concretely modeling the effect of microplankton dynamics on mesozooplankton grazers. Addressing them systematically in models will require novel integration between what could be called oceanographic and marine-biological perspectives on large zooplankton.

CONCLUSION
Coltrane 1.0, introduced here, is a minimalist model of copepod life history and population dynamics, a metacommunity-level framework on which additional species-or population-level constraints can be layered. Many present and future patterns in large copepods might well prove to be sensitive to species-specific constraints that Coltrane 1.0 does not resolve, such as thermal adaptation, physiological requirements for egg production, or cues for diapause entry and exit. Nevertheless, the model experiments above demonstrate that many patterns in latitude, time, and trait space can be replicated numerically even when we only consider a few key constraints on the individual energy budget: the total energy available in a given environment per year; the energy and time required to build an adult body; the metabolic and predation penalties for taking too long to reproduce; and the size and temperature dependence of the vital rates involved.
Results of the global and Bering experiments (Figures 5-7) suggest that timing and seasonality are crucial to large copepods, but not because of match/mismatch (Edwards and Richardson, 2004) the model organisms are free to resolve timing mismatches with complete plasticity. Rather, these results highlight the role of seasonality in the sense of total energy available for growth and development per year, or the number of weeks per year of net energy gain relative to the number of weeks of net deficit. The simplicity of this view means that the model scheme and results may generalize far beyond copepods with only minor modification.
The exercise of parameterizing the Bering Sea and Disko Bay cases, and of attempting to map real environments onto an idealized parameter space in the global experiment (Figure 6), highlighted that the real limit on our ability to predict the fate of copepods in changing oceans may not be our incomplete knowledge of their physiology, but rather our incomplete knowledge of how their environments appear from their point of view. How do standard oceanographic measures of chlorophyll and particulate chemistry relate to prey quality, and how much risk a copepod should take on in order to forage in the euphotic zone? How do bathymetry, the light field, and other metrics relate to the predator regime? Further experiments in a simple, fast, mechanistically transparent model like Coltrane may suggest new priorities for field observations, in addition to new approaches to regional and global modeling.

AUTHOR CONTRIBUTIONS
NB designed the model, performed the analysis, and led the writing of the manuscript. EM and TN helped formulate and interpret the Disko Bay case study, and LE the Bering Sea case study. All authors contributed to revision of the manuscript.

FUNDING
This work was supported by grants PLR-1417365 and PLR-1417224 from the National Science Foundation.