Journal of Plankton Research

JPR Advance Access originally published online on June 21, 2009
Journal of Plankton Research 2009 31(9):939-963; doi:10.1093/plankt/fbp042

This Article Abstract FREE Full Text (PDF) Supplementary Data All Versions of this Article: 31/9/939    most recent fbp042v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Google Scholar Articles by Cropp, R. Articles by Norbury, J. PubMed Articles by Cropp, R. Articles by Norbury, J. Social Bookmarking          What's this?

© The Author 2009. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oxfordjournals.org

Parameterizing plankton functional type models: insights from a dynamical systems perspective

Roger Cropp^1,* and John Norbury²

¹ Centre for Environmental Systems Research, Griffith School of Environment, Griffith University, Nathan, QLD 4111, Australia ² Mathematical Institute, University of Oxford, 24-29 St Giles, Oxford OX1 3LB, UK

^* CORRESPONDING AUTHOR: r.cropp{at}griffith.edu.au

Received on March 31, 2009; accepted on May 22, 2009

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
The spectre of anthropogenic global climate change has focused attention on biogeochemical cycling in the oceans as marine plankton ecosystems are involved in the cycling of several compounds thought to have significant implications for climate. To better understand these processes, modellers are developing plankton functional type (PFT) models that group plankton according to their biogeochemical properties. There is some debate as to whether our understanding of plankton ecosystems is sufficiently well developed for PFT models to be reliable and for their predictions to be treated with confidence. In this paper, we examine the dynamical properties of a generic predator–prey–prey PFT model, then apply these analysis techniques to a simple example PFT model with two phytoplankton and one zooplankton in order to explore its parameter space. We find that parameter combinations for which all PFTs stay extant for all time appear rare, but develop a simple heuristic that allows such parameter sets to be identified relatively easily for many PFT models. We observe that such systems often have phytoplankton with similar growth rates, but that differ in other properties such as differing nutrient utilization strategies or different susceptibilities to grazing. We also note that persistent PFT systems are more likely if neither phytoplankton have a low specific mortality rate or is a highly nutritious food for the grazer.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
The prospect of anthropogenic climate change has generated substantial interest in the role of plankton in biogeochemical cycling in the ocean. Plankton may have a significant influence on climate by drawing down carbon dioxide from the atmosphere, sequestering it in the deep ocean, and by producing dimethylsulphide and other volatile compounds that may affect cloud formation over the oceans. Plankton models that include several plankton functional types (PFTs) are needed to resolve the role of plankton in biogeochemical cycling, as different plankton utilize different elements in different ways. PFTs are therefore often defined according to their biogeochemical role rather than their phylogeny (Falkowski et al., 2003). The number of PFTs required is dependent on the purpose of the model, with Hood et al. (Hood et al., 2006) noting that the performance of PFT models is more closely related to their tuning than their complexity. Le Quéré et al. (Le Quéré et al., 2005) suggested that 10 key PFTs were required to resolve climatically important biogeochemical cycling in the oceans. However, there has been significant debate over whether plankton ecosystems are sufficiently well understood to place any reliance on the results of models that include multiple PFTs (Anderson, 2005; Le Quéré, 2006).

Pragmatically, the demand to resolve the roles of different PFTs in biogeochemical cycling in the oceans has meant that PFT models are already being developed and applied. The debate over the merits of including phytoplankton functional types in the generally successful nutrient–phytoplankton–zooplankton (NPZ) models suggests a pressing need for better understanding of the behaviour of these models. The emphasis of attempts to improve the understanding of, and build confidence in, PFT models is often focussed on more and more accurate measurements of PFT traits (parameter values) and more and more accurate data to calibrate and validate the models (Le Quéré et al., 2005; Hood et al., 2006). Here, we investigate the dynamics of a simple, generic model with two phytoplankton functional types to analyse what insights might be gained into the attributes of more complex PFT models. We choose the model and its parameters to ensure ecological realism under the conditions derived by Kolmogorov (Kolmogorov, 1936) and explicated by May (May, 1973).

Many climate change computer simulations include a model of plankton ecosystem processes. Many of these models of plankton dynamics, both simple and complex (Spitz et al., 2001; Franks, 2002; Cropp et al., 2004; Vallina et al., 2008), may be classed as Kolmogorov systems, as they are of the general form:

(1)

where for t > 0, and the functions f_i are bounded and continuously differentiable in their variables u₁, u₂, ... , u_n. The f_i describe the growth and mortality of each species, trophic guild or PFT, that is f_i = (growth – predation – mortality)_i. These f_i are often nonlinear (but smooth) functions of u₁, u₂, ... , u_n and include several parameters that describe the attributes (traits) of the plankton and how they interact. There are many options for the f_i, and much of the debate surrounding the application of PFT models centres on whether the forms of the f_i and the values of the parameters that distinguish one PFT from another are sufficiently well understood. However, the f_i commonly used in ecosystems models typically meet the above requirements.

In addition to defining these systems, Kolmogorov (Kolmogorov, 1936) derived conditions on the f_i for predator–prey interactions in the form of equation (1) for n = 2 that ensured only ecologically realistic dynamics (stable spiral equilibriums or stable limit cycles that ensure continued co-existence of both predator and prey) were possible. Systems with stable spiral equilibrium points will exhibit blooms of ever-decreasing amplitude until they achieve a steady state, whereas systems with stable limit cycles will eventually settle to periodic repeating blooms of identical amplitude.

Ecological interpretations of Kolmogorov's conditions are provided by May (May, 1973). A consequence of Kolmogorov's conditions is that they require the system to have both an autotroph-only and a predator–prey critical (equilibrium) point in the ecologically feasible region of the state space where all state variables are positive. The existence of these critical points is determined by the nature of the f_i, but whether they lie within the ecologically feasible region of the state space is determined by the parameter values used in the model.

We consider the endogenous dynamics of a model with multiple phytoplankton functional types that is a three-dimensional (n = 3) Kolmogorov system.

(2)

where u₁ and u₂ represent phytoplankton functional types and u₃ represents their zooplankton grazer. As the f_i are generally nonlinear functions of u₁, u₂, and u₃, analytic solutions to such systems are rare. We shall restrict our analysis to f_i that comply with Kolmogorov's (Kolmogorov, 1936) conditions, as do many ecosystem models in the contemporary literature (Huang and Zhu, 2005).

We look at the particular case of a model that conserves the mass of limiting nutrient as many models applied in biological oceanography also have this property (Spitz et al., 2001; Franks, 2002; Vallina et al., 2008). Conservation of mass implies that the total mass of inorganic nutrient (N) present at any time is given by:

(3)

where N_T is a constant that gives the total effective nutrient in the system and the u_i are the concentrations of the PFTs measured in this currency. We note that conservation of mass is required for a biogeochemical model to be written as a Kolmogorov system, as the nutrient equation in these models typically cannot be written in the Kolmogorov form. However, conservation of mass causes N to become a "virtual" variable; as shown in equation (3), it allows the equation to be inferred from the other equations that are in Kolmogorov form.

We scale each state variable with respect to the total nutrient (i.e. u_i/N_T), which allows us to define an ecologically feasible "state space" where 0 u₁, u₂, u₃ 1. Conservation of mass then implies that the dynamics of the system are confined to the part of a multi-dimensional Cartesian co-ordinate space where each axis represents the concentration of one species or functional type (such as u₁, u₂, u₃) and all the variables are positive and <1. Mass conservation is said to apply when on N = 0, i.e. this gives the crucial physically realistic constraint on our f_i, which then forces the dynamics to live in the feasible state space with u₁ + u₂ + u₃ 1 for all time.

We consider the critical (equilibrium) points of this system, denoted by where for all time. Implicit in the rationale for constructing plankton models with more than one functional type is the assumption that, in the absence of environmental factors, an interior critical point, with , both exists and is an important determinant of the dynamics of the system. The theme of this work is an inquiry into the nature of these interior (predator–prey–prey) critical points in systems with two autotrophs and a grazer, and how this might inform the development and calibration of more complex PFT models. We will show that most parameter sets do not produce useful PFT models; carefully tuned models are essential. While the range of parameter values may be narrowed down by laboratory studies, as Le Quéré (Le Quéré, 2006) observes, our results suggest that laboratory studies need to include estimates of parameters for combinations of organisms, as suggested by Flynn (Flynn, 2006). We also note that, while parameter sets that produce realistic PFT models are rare, they are distributed throughout the parameter space, suggesting that in addition to laboratory measurements, model dynamics should also be rigorously validated against the observed data (Anderson, 2006). We derive biological heuristics for parameter combinations that assist in the tuning of parameter sets and yield useful PFT models.

	METHOD

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
We first analyse the dynamics of a generic three state variable Kolmogorov system (equation (2)) where the u₁ and u₂ represent autotrophs and u₃ represents their grazer. The analysis of this generic system provides general results that apply to all such three state variable Kolmogorov systems, irrespective of the process formulations (f_i) chosen to represent the interactions between the state variables and independent of the parameter values used in the model.

We consider the critical points of the system, defined by for all i. In Kolmogorov systems, critical points may be obtained from the isoclines in two ways for each equation, when f_i = 0 or when u_i = 0. Each critical point in a three state variable system has three eigenvalues and three associated eigenvectors. These eigenvalues and eigenvectors describe the local (Lyapunov) stability of the system in the region of each of the critical points, and together form "signposts" that control the dynamics of the system. Negative eigenvalues attract nearby trajectories (states of the system) to the critical point along their eigenvectors, whereas positive eigenvalues repel nearby trajectories. Critical points that have all of their eigenvalues negative are stable, and produce systems that are attracted to, and then remain in, equilibrium. However, it only requires one eigenvalue of a critical point to be positive to render the point unstable, and the system will never be in the state described by the critical point for more than a brief period. The properties of critical points determined by their eigenvalues and eigenvectors will be discussed in more detail below.

The isoclines, critical points and eigenvectors for the whole system are shown in Fig. 1. Some eigenvectors have been omitted for clarity; none are shown for point E but these can be inferred from those shown for D and F. Eigenvectors that always have positive eigenvalues are shown pointing away from the critical point, those that always have negative eigenvalues are shown pointing towards the critical point, and the others have double-ended arrows. Figure 1a shows the three-dimensional state space in which the isoclines are surfaces; where parts of the isoclines are hidden behind other isoclines they are shown by dotted lines. Figure 1b shows the predator–prey state spaces ( and ) that form the left and right faces of the three-dimensional space in Fig. 1a, while Fig. 1c shows the competitive autotroph state space that forms the base of Fig. 1a. The dashed lines indicate the conservation of mass conditions for each face.

View larger version (24K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 1. Generic diagram of the zero isoclines (f1, f2, f3 = 0 and u1, u2, u3 = 0), critical points (stars) and eigenvectors (bold arrow) for the Kolmogorov system described by equation (2). (a) The full system (u1, u2, u3), (b) the predator–prey subsystems (u1, u3 and u2, u3) and (c) the competitive autotroph subsystem (u1, u2).

 
The eigenvalues and eigenvectors of a system reflect the essential system-level characteristics of the community, or Jacobian, matrix that describes how each variable changes in response to changes in itself and other variables. Technically, the elements of the community matrix (J) are given by the partial derivatives of the equations with respect to the state variables, that is J_i,j = (f_i u_i)/u_j for i, j = 1,2, ... , n.

The critical points and their associated eigenvalues and eigenvectors control the dynamics of the system by influencing the trajectories that the system takes through the state space. The state space of a three state variable Kolmogorov system is a three-dimensional Cartesian coordinate space with u₁, u₂ and u₃ along the axes (Fig. 1). A trajectory is defined by an initial condition (i.e. the initial values ascribed to each of u₁, u₂ and u₃) and the subsequent behaviour of the system over time as defined by the equations . A negative eigenvalue at a critical point means that trajectories near the eigenvector associated with the eigenvalue are attracted to the critical point along the eigenvector. Conversely, a positive eigenvector at a critical point indicates that trajectories near the associated eigenvector will be repelled from the critical point along the direction of the eigenvector. Each critical point has as many eigenvalues as the dimension of the system (in this case three). Critical points that have all negative real parts are locally stable, that is, all nearby trajectories will be attracted towards them and the system will display an equilibrium state, but it only requires one positive eigenvalue for a point to be unstable.

Eigenvalues that are real numbers indicate that nearby trajectories will approach/retreat from the critical point asymptotically (i.e. continuously increasing or continuously decreasing). However, many predator–prey systems have eigenvalues that are complex numbers, with a real part and an imaginary part. The real part of the eigenvalue determines whether nearby trajectories are attracted to, or repelled from, the critical point, while the complex part indicates that nearby trajectories will spiral around the critical point (i.e. the system will oscillate). We shall use the following terms to describe the stability properties of the critical points:

• Stable: the point has three negative real eigenvalues
  
• Unstable or saddle: the point has at least one positive real eigenvalue
  
• Stable spiral: the point has a pair of complex eigenvalues that have negative real parts
  
• Unstable spiral: the point has a pair of complex eigenvalues that have positive real parts.
  

When we discuss the stability properties of the predator–prey–prey critical point, which typically has a pair of complex eigenvalues and a single real eigenvalue, we shall differentiate its stability properties into those of the complex eigenvalues and those of the real eigenvalue.

We then consider a specific example of an NP₁P₂Z Kolmogorov system in order to examine the parameterizations of the system that result in realistic PFT dynamics. The study system has conventional phytoplankton (P₁, P₂) growth on inorganic nutrient (N) balanced by zooplankton (Z) grazing and linear mortality. We note that there is currently significant debate over the appropriate form of the process representations in PFT models (Flynn, 2003; Mitra, 2009). We choose "simple" f_i that are commonly used to allow, as far as possible, for closed form analytic expressions to be found for the critical points and their eigenvalues. Even in this simple model, explicit analytic evaluation is not always possible, and we are forced to develop numerical estimates of one interior value. We justify this choice because

• we will demonstrate in the analysis of the generic Kolmogorov system that many of our results are independent of the form of the f_i, and
  
• at the heart of our approach is the great value that analytic expressions provide over numerical estimates in understanding why PFT models behave the way they do.
  

Analytic expressions for the key properties of PFT systems allow us to understand and predict PFT dynamics rather than merely observe them. We shall present the results of the analysis of the NP₁P₂Z model, equivalent to our analysis of the generic Kolmogorov system, and use these results to explore the parameter space and associated dynamical properties of this simple system.

The example system is written in a currency of inorganic nutrient, with all state variables expressed as concentrations of nutrient, as described by equations (4)–(6):

(4)

(5)

(6)

We check our conservation of mass criterion as per equation (3) and see that

(7)

when we put N = 0 in the right-hand side then at N 0 = N_T–P₁– P₂– Z. This forces N to be positive, so that P₁ + P₂ + Z < N_T for all time in our system. The useful PFT dynamics, where nothing goes extinct, occur in 0 < P₁, P₂, Z with P₁ + P₂ + Z < N_T. We shall define several parameter sets for this model to demonstrate qualitatively different dynamics. These parameter sets are derived from measured values that are described in Table I. As we will be considering the effect of varying parameter values on the model dynamics, it is essential that we scale the parameters so that the influence of each parameter is revealed unequivocally. We have therefore non-dimensionalized the parameter values in Table I from their measured values by scaling time by the maximum growth rate of P₁(µ) and concentrations by the total nutrient (N_T).

View this table: [in this window] [in a new window]   Table I: Measured parameter values for the NP1P2Z model

 
We have chosen the form of the Z grazing terms so that the parameter values of and in equations (4)–(7) can be used to control the dynamics of the predator–prey sub-systems (NP₁Z and NP₂Z). Using values of or equal to zero results in these terms behaving as Lotka–Volterra grazing formulations, while non-zero values cause them to behave as Michalis–Menten grazing formulations. This provides a simple mechanism to determine the dynamics of the predator–prey subsystems as it is easy to find parameter sets that result in asymptotically stable spiral dynamics in this system with Lotka–Volterra grazing, and similarly easy to find parameter sets that result in stable limit cycle dynamics in this system with Michalis–Menten grazing. In general, in what follows we will set 0 0.1 to generate stable spirals in the NP₁Z subsystem and to generate stable limit cycles in the NP₂Z subsystem. This system can therefore reproduce all of the dynamical states that Kolmogorov (Kolmogorov, 1936), May (May, 1973) and most subsequent ecologists have considered ecologically realistic for predator–prey systems.

We consider the dynamics of the NP₁P₂Z system under six parameter sets that provide six representative dynamical behaviours, and examine the influence of the various eigenvalues and eigenvectors on the dynamics. In each case, typical trajectories are integrated from two different initial conditions with a fourth–fifth-order adaptive step size Runge–Kutta routine implemented with relative and absolute step size tolerances set to machine epsilon. The dynamics of the system for the six parameter sets explicate the dynamics of three classes of PFT systems: bona fide PFT systems where all functional types remain extant forever; pseudo-PFT systems where the initial population sizes determine which type will go extinct; and non-PFT systems where the parameter values determine which type will go extinct irrespective of the initial conditions.

We then investigate the sensitivity of the interior predator–prey–prey critical point of a bona fide PFT system to variations in parameter values. This analysis is not intended to be a comprehensive sensitivity analysis of the model, but a demonstration of the finely balanced nature of bona fide PFT systems, and an explication of the important parameter sensitivities of the example model at one point in the parameter space (column 1 of Table II). The parameters were varied one at a time and the location of the predator–prey–prey critical point was followed using a numerical continuation algorithm (Dhooge et al., 2003), which follows the critical point through the state space, until either the parameter or the critical point achieved ecologically infeasible values. The maximum and minimum values of the parameter that resulted in the predator–prey–prey critical point being in the feasible region of the state space were recorded, and a sensitivity metric calculated by scaling the range of valid parameter values by the mean parameter value. The parameter set used resulted in a stable spiral on one face of the state space, characterized by eigenvalues with negative real parts, and a limit cycle on the other, characterized by eigenvalues with positive real parts. We therefore also recorded the parameter value at which the real parts of the eigenvalues passed through zero and the dynamics changed from a stable spiral to a limit cycle or vice versa as the predator–prey–prey critical point traversed the state space.

View this table: [in this window] [in a new window]   Table II: Parameter values used in the numerical simulations

 
Finally, we consider the ubiquity of parameter sets that result in bona fide PFT systems. The biological parameter spaces that marine biogeochemical models occupy are notoriously poorly constrained, so we defined a parameter space ranging from 50% of the smallest to 150% of the largest values of each parameter listed in Table II. We randomly sampled 5 x 10⁶ parameter sets from uniform distributions within this parameter space. We established criteria for the validity of each parameter set based on the application of Kolmogorov's (Kolmogorov, 1936) criteria to the analysis of our exemplar model. If a parameter set was valid, we calculated the eigenvalues of the inward-pointing eigenvectors of the predator–prey critical points on the faces and classified each parameter set according to the signs of the eigenvalues. We tested the effect of extending the parameter space beyond the limits specified above, but this did not significantly alter the number or proportion of valid parameter sets.

View larger version (40K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 2. Critical points and dynamics for the NP1P2Z system showing attraction of both trajectories to an internal stable spiral critical point. Parameter set used is listed in column 1 of Table II. Initial conditions for red trajectory are N = 0.898, P1 = 0.100, P2 = 0.001, Z = 0.001. Initial conditions for green trajectory are N = 0.799, P1 = 0.100, P2 = 0.100, Z = 0.001.

 

View larger version (42K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 3. Critical points and dynamics for the NP1P2Z system showing attraction of both trajectories to an internal stable limit cycle. Parameter set used is listed in column 2 of Table II. Initial conditions for red trajectory are N = 0.898, P1 = 0.100, P2 = 0.001, Z = 0.001. Initial conditions for green trajectory are N = 0.655, P1 = 0.092, P2 = 0.092, Z = 0.161.

 

View larger version (38K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 4. Critical points and dynamics for the NP1P2Z system showing attraction of both trajectories to a stable spiral on each face. Parameter set used is listed in column 3 of Table II. Initial conditions for red trajectory are N = 0.849, P1 = 0.100, P2 = 0.050, Z=0.001. Initial conditions for green trajectory are N = 0.849, P1 = 0.050, P2 = 0.100, Z = 0.001.

 

View larger version (39K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 5. Critical points and dynamics for the NP1P2Z system showing attraction of the red trajectory to a stable spiral on the P1 face and the green trajectory to a limit cycle on the P2 face. Parameter set used is listed in column 4 of Table II. Initial conditions for red trajectory are N = 0.899, P1 = 0.050, P2 = 0.050, Z = 0.001. Initial conditions for green trajectory are N = 0.849, P1 = 0.050, P2 = 0.100, Z = 0.001.

 

View larger version (57K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 6. Critical points and dynamics for the NP1P2Z system showing attraction of both trajectories to the stable spiral on the P1 face. Parameter set used is listed in column 5 of Table II. Initial conditions for red trajectory are N = 0.799, P1 = 0.100, P2 = 0.100, Z = 0.001. Initial conditions for green trajectory are N = 0.889, P1 = 0.010, P2 = 0.100, Z = 0.001.

 

View larger version (53K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 7. Critical points and dynamics for the NP1P2Z system showing attraction of both trajectories to the stable limit cycle on the P2 face. Parameter set used is listed in column 6 of Table II. Initial conditions for red trajectory are N = 0.799, P1 = 0.100, P2 = 0.100, Z = 0.001. Initial conditions for green trajectory are N = 0.889, P1 = 0.100, P2 = 0.010, Z = 0.001.

 

	RESULTS

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
We initially present the results of the analysis of the generic Kolmogorov system (equation (2)), where we define u₁ and u₂ to be autotrophs and u₃ to be their grazer. We present analytic expressions for the critical points and eigenvalues where possible, and present the generic dynamics of the state space described by the critical points, eigenvalues and eigenvectors graphically (Fig. 1). We have labelled the critical points in Fig. 1 and shall use these labels to identify the critical points and their eigenvalues in the following analysis. The values of the state variables at the critical points are denoted by the asterisk superscripts, and the critical points differentiated by their subscripts. The signs of the eigenvalues are shown where these are defined by the ecological properties of the system and always hold.

Origin critical point (O)

Every Kolmogorov system has a critical point at the origin (O in Fig. 1) where (i.e. no life) and (usually) f_i 0 for all i. The eigenvalues at the origin are

(8)

(9)

(10)

where means that the expression f_i is evaluated at the critical point O (i.e. using the values of the state variables at the critical point). Analytic expressions for the eigenvalues that control the stability characteristics of this point are therefore easily obtained from inspection of the equations as they are simply the f_i of each PFT, and their magnitudes given by the value of the f_i at the origin.

The origin represents the state of the system where only inorganic nutrient exists. Near this point, autotrophs must grow (or there could never be any life) and predators must die (as they have insufficient prey to survive on). This point will therefore always be a saddle (i.e. have some positive and some negative eigenvalues) and its eigenvectors (shown by the bold arrows in the Fig. 1) will always lie along the axes of the state space. The unstable directions (_{O –1} and _{O –2}) are a result of the autotrophs growing by consuming nutrient in isolation along their respective axes, while the stable direction (_{O –3}) is that of the predator dying in the absence of prey along its axis. Note that each eigenvalue is comprised only of parameters related to each PFT's own growth and/or mortality, and that grazing interactions are not relevant to this critical point. Further, ecological realism requires that it should be possible to make the origin a stable node. The f_i of the autotrophs should therefore all include an explicit mortality term such that it is possible to increase the mortality rates of the autotrophs so that they exceed the maximum growth rates possible and extinction of all life ensues.

Prey-only critical points (A, C)

Every model that complies with Kolmogorov's (Kolmogorov, 1936) criteria will have a prey-only (autotroph) critical point in each of the predator–prey subsystems. These are defined by (A) and (C) and are located where the f₁ = 0 isocline in the (u₁, u₃) plane intersects the u₁ axis and where the f₂ = 0 isocline in the (u₂, u₃) plane intersects the u₂ axis, respectively. The eigenvalues of these points are given by, for A:

(11)

(12)

(13)

and for C:

(14)

(15)

(16)

Here again the eigenvalues have neat ecological interpretations. The stable eigenvalues (_{A –1} for A and _{C –2} for C) are given by the response of the blooming autotroph to increases in its own biomass (f_i/u_i) and are always negative, as autotroph growth rates reduce as nutrient becomes less available. The eigenvectors associated with these eigenvalues always lie along the u₁ and u₂ axes, respectively, as shown by the bold arrows in Fig. 1.

The other eigenvalues are obtained by evaluating the growth functions (f_i) of the other PFTs at the critical point. These reflect the grazing pressure applied by the predator (_{A –3} for A and _{C –3} for C) or the competition from the other autotroph for nutrient (_{A –2} for A and _{C –1} for C). Systems that comply with Kolomogorov's criteria will always have _{A –3} and _{C –3}, the eigenvalues associated with grazing pressure, positive (destabilizing) at these points. The directions of these eigenvalues will vary according to the nature of the f_i, but will always point into the interior of, and lie in the plane of, the (u₁, u₃) plane for A or the (u₂, u₃) plane for C. The eigenvalues associated with the competing autotroph (_{A –2} for A and _{C –3} for C) may be positive if the blooming autotroph leaves enough nutrient so that the competing autotroph's growth from nutrient uptake exceeds its specific mortality rate, or negative otherwise. The eigenvectors for these critical points are always orthogonal to the (u₁, u₃) plane for A or the (u₂, u₃) plane for C as shown in Fig. 1.

Dual prey-only critical point (B)

Figure 1 shows a third prey-only critical point (B), where both competing autotrophs coexist. The existence of this point is not assured by Kolomogorov's criteria, as the (u₁, u₂) plane does not constitute a predator–prey system. This point is defined by and is located where the f₁ = 0 isocline intersects the f₂ = 0 isocline in the (u₁, u₂) plane. The eigenvalues of this point are given by:

(17)

(18)

The eigenvalue associated with grazing pressure (_{B –3}) is always positive, whereas the eigenvalues associated with the autotroph growth (_{B –1,2}) will have one negative eigenvalue representing autotroph growth on available nutrient. The other eigenvalue may be positive or negative, either repelling or attracting nearby trajectories, but is not critical to the co-existence of u₁ and u₂. The directions of the eigenvectors associated with these eigenvalues will vary according to the f_i and the parameter values used but will generally lie as shown in Fig. 1 for parameter values that result in B being located in the feasible region of the state space.

Predator–prey critical points (D, F)

Every sub-system that complies with Kolmogorov's criteria will have a predator–prey critical point. In the system defined by equation (2) these are defined by for D and for F and are located where the f₁ = 0 isocline intersects the f₃ = 0 isocline in the (u₁, u₃) plane and where the f₂ = 0 isocline intersects the f₃ = 0 isocline in the (u₂, u₃) plane, respectively. The eigenvalues of these critical points are given by:

(19)

(20)

for D, and by:

(21)

(22)

for F. In almost all cases _{D –1,3} and _{F –2,3} will be complex numbers, with positive or negative real parts, indicating that trajectories will either spiral into or away from the critical point. An example of each case is shown in Fig. 1, with D shown with negative real parts and F shown with positive real parts. Spiral curves then lie in the (u₁, u₃) or (u₂, u₃) planes and start or end at D or F. _{D – 1,3} and _{F – 2,3} therefore control the dynamics in the (u₁, u₃) and (u₂, u₃) planes, respectively.

The dynamics of the system in the direction orthogonal to these planes is of critical importance to PFT modellers as the eigenvalues in this direction determine whether a system will maintain all plankton extant during simulations. The eigenvalues in this direction are given by _{D –2} for D and by _{F –1} for F. The eigenvalues are associated with eigenvectors that are orthogonal to the (u₁, u₃) and (u₂, u₃) planes, respectively, with magnitudes given by the f_i of the competing autotroph evaluated at the critical point. These eigenvalues are always real numbers, their direction is known, and analytic expressions for them are easily obtained from inspection of the model equations. The signs of these eigenvalues determine whether a predator–prey–prey critical point (E), fundamental to the construction of a bona fide PFT model, exists and is stable.

If _{D –2} and _{F –1} are both positive, E exists in the feasible region of state space and is stable in the direction orthogonal to the predator–prey planes (u₁, u₃) and (u₂, u₃). This means that the system will eventually come to some balance where u₁ and u₂ co-exist, and are both grazed on by u₃. This may be a stable equilibrium point or a stable limit cycle.

If _{D –2} and _{F –1} are both negative, E exists in the feasible region of state space but is unstable in the direction orthogonal to the predator–prey planes (u₁, u₃) and (u₂, u₃). This means that u₁ and u₂ cannot co-exist, and one must always out-compete the other. The winner and loser of the competition in this case are determined by their starting values (initial conditions). We shall refer to these cases as pseudo-PFT systems. When _{D –2} and _{F –1} have opposite signs, E does not exist in the feasible region of state space, and again u₁ and u₂ cannot co-exist. In this case, the winner and loser are determined by the eigenvalues, and the initial conditions have no influence on the outcome of the competition. We shall refer to these cases as non-PFT systems.

The ecological interpretation of either _{D –2} or _{F –1} being positive is that if the equilibrium state (or equivalently the long term average if a limit cycle) of the competing autotroph and its grazer leaves enough nutrient available for the autotroph growth rate to exceed its losses due to grazing and mortality, then the autotroph has the capacity to grow. Technically, this destabilizes the critical point allowing the trajectory to explore the interior of the state space.

Predator–prey–prey critical point (E)

As noted above, the system may have a critical point E defined by located where the f₁ = 0, f₂ = 0 and f₃ = 0 isoclines all intersect in the (u₁, u₂, u₃) volume. In this case the isoclines are surfaces, as shown in Fig. 1. The existence of this point is not assured by Kolmogorov's criteria as it is a predator–prey–prey system rather than a predator–prey system. However, systems that have critical points at D and F will, for some parameter sets, have a predator–prey–prey critical point E that lies in the interior of the state space.

The eigenvalues of the critical point E are generally intricate in analytic form, and difficult to interpret, as they involve the roots of a cubic equation derived from the community matrix, and hence it is usually simpler to obtain them numerically. It appears that the spiral dynamics enforced on the (u₁, u₃) and (u₂, u₃) planes by the Kolmogorov criteria are reflected throughout the interior of the state space. We therefore expect that the critical point E will have one pair of complex conjugate eigenvalues that control its spiral behaviour in the (₁P₁ + ₂P₂, Z) "plane", where ₁ and ₂ are constants defining a ray on the (P₁, P₂) face. The third eigenvalue must be a real number, and its sign will control the outcome of the competition between the two autotrophs. The real eigenvalue will therefore either repel trajectories away from E toward the predator–prey critical points on the faces (if positive) or attract them away from the faces toward the interior predator–prey–prey critical point E (if negative).

The global dynamics of the state space must be consistent with the above local information. Therefore, if the interior critical point E has positive real eigenvalues, the critical points on the faces D and F must both have negative real eigenvalues associated with the eigenvectors that are orthogonal to the faces. Similarly, if E has negative real eigenvalues, D and F must both have positive real eigenvalues. In the case where D and F have real eigenvalues of opposite signs, E cannot exist in the interior of the state space. In this circumstance, the point is located exterior to the ecologically feasible state space, and it and its eigenvalues do not influence the dynamics of the system. These cases will be demonstrated in our analysis of the example NP₁P₂Z Kolmogorov system to follow.

Critical points of the NP₁P₂Z system

Although relatively simple, the NP₁P₂Z system lies at the boundary of complexity for systems that are wholly analytically tractable. We are able to obtain analytic expressions for all but one of the critical points and eigenvalues; the location of the interior predator–prey–prey critical point E and its eigenvalues and eigenvectors is easily obtained numerically. The analytic expressions describing the locations of the critical points, and their eigenvalues where possible, analogous to and labelled identically to that presented above for the generic model are presented in the Appendix. The numerical evaluations of these expressions are presented in Table III, and the locations of the critical points that these values specify are shown in Figs 2–7. Note that critical point B is a special case of a critical line in this model (see Appendix for details) and is therefore not shown in Figs 2–7 that show the dynamics of the system.

View this table: [in this window] [in a new window]   Table III: Summary of critical point locations and their eigenvalues for the parameter values used in the numerical simulations

 
We now review the dynamics possible in NP₁P₂Z systems. This is not a comprehensive review, but rather a pragmatic survey of the dynamics; other dynamics can be inferred from the sample presented. We formulated six parameter sets (Table II) derived from the measured values listed in Table I. These parameter sets were chosen to produce dynamics representative of the six qualitatively different cases, and always result in a complex conjugate pair of eigenvalues in the (P₁, Z) or (P₂, Z) vertical planes, whereas the eigenvalues in the (P₁, P₂) horizontal plane are always real. The dynamics shown in Figs 2–7 demonstrate the various possibilities of the signs of the real parts of these eigenvalues. The figures each show two trajectories that start from different initial conditions to demonstrate the different trajectories possible. The arrows on the faces represent the vector fields of the sub-systems (i.e. NP₁Z, NP₂Z and NP₁P₂). The vector fields are determined by the eigenvalues and eigenvectors on the faces and allow generalized dynamics to be inferred from the individual trajectories shown. Colour versions of these figures are provided online in the Supplementary material.

Parameter sets for which both real eigenvalues of the predator–prey critical points on the faces are positive are the only cases in which the original set of functional types is always maintained in the model. In all other cases (i.e. both negative or one positive and one negative eigenvalue), one species or other of the phytoplankton will dominate. The system dynamics will then no longer inhabit the interior of the state space, but will reside on either the (P₁, Z) or (P₂, Z) face.

We observe that linear constant coefficient ordinary differential equation models have exponential solutions. For ecological models, this means that, in the absence of growth, state variables with linear mortality terms exponentially decay. They cannot therefore ever become exactly equal to zero (although when solved numerically they may equal zero, or even become negative, if the numerical integration routine used is not sufficiently accurate). However, we note that when the dynamics of our PFT system are attracted to one of the faces then the other PFT, even though non-zero, can never recover, and is for all biogeochemical purposes extinct. For ease of discussion, we will therefore refer to this situation as the extinction of that PFT.

Bona fide PFT systems

We define bona fide PFT systems to be those in which no PFT goes extinct. There is therefore a consistency between the model equations that describe multiple PFTs and the numerical simulations that maintain those PFTs for all time. Figure 2 presents typical dynamics for a bona fide PFT system with a stable spiral internal critical point. The origin and autotroph critical points are saddles and the predator–prey critical points on the faces both have positive real eigenvalues (equations (A.26) and (A.27)). The complex eigenvalues of the predator–prey critical point on the (P₁, Z) face have negative real parts, whereas the complex eigenvalues of the predator–prey critical point on the (P₂, Z) face have positive real parts. The (P₁, Z) face therefore has stable spiral dynamics, whereas the (P₂, Z) face has stable limit cycle dynamics.

Every initial condition for the parameterization of the model in Fig. 2 will end up at the spirally stable predator–prey–prey critical point in the interior of the feasible state space. This is perhaps the ideal dynamics for a PFT model, as a balance between all PFT is achieved, and no PFT is ever lost to the system.

Figure 3 demonstrates the dynamics of perhaps the next most desirable system for PFT modelling, that of a stable limit cycle predator–prey–prey critical point. Again, the (P₁, Z) face has stable spiral dynamics while the (P₂, Z) face has stable limit cycle dynamics, and both have positive real eigenvalues. Here, the system regularly cycles on a plane in the interior of the state space, in this case with large blooms of P₂ and relatively small blooms of P₁. Again, these dynamics are robust and no PFT is lost to the system as the real eigenvalues of the predator–prey critical points on the faces direct all orbits away from the faces and towards the interior point.

Pseudo-PFT systems

We define pseudo-PFT systems to be one class of models for which there is an inconsistency between the model that is ostensibly described by the system of differential equations and the results of the numerical solution of the model. In our pseudo-PFT systems, one or other of the phytoplankton functional types will always go extinct in simulations, with the PFT destined for extinction being determined by the initial conditions selected to commence the integration from. Figure 4 shows the first of the cases in which extinction of one P occurs. Both (P₁, Z) and (P₂, Z) faces have stable spiral dynamics in this case, and also have negative real eigenvalues attracting nearby orbits "horizontally" onto the face. This is perhaps the worst scenario for the dynamics of a PFT model, as the initial dynamics may give the impression of a bona fide PFT system. However, the system will always collapse to one or the other of the stable states on the faces, and which P survives is determined solely by the initial conditions. We note, however, that for different parameter sets that result in this dynamical situation, the survivorship of the P will be determined by different initial conditions, that is, the basins of attraction of the predator–prey critical points on the faces will be different for different parameter sets. We further observe that it would be possible to start simulations with phytoplankton functional type populations very close to the separating surface dividing dominance by P₁ from dominance by P₂. In such an event, the system would give the initial impression that it was maintaining both types, but long simulations would eventually result in the extinction of one functional type.

Figure 5 shows a similar situation to Fig. 4, as both predator–prey critical points have negative real eigenvalues, except in this case the (P₂, Z) face has stable limit cycle dynamics. It is interesting to note that although the initial conditions again solely determine which PFT will survive, it is not necessarily immediately obvious which PFT that will be. The left trajectory, for example, reveals that P₁ dominates the initial bloom before the trajectory sweeps over to near the (P₂, Z) face, where subsequent smaller blooms are dominated by P₂. As subsequent blooms wax and wane, however, P₁ gradually re-asserts its dominance, resulting in the "extinction" of P₂ in this case.

Non-PFT systems

We define non-PFT systems to be those in which one phytoplankton functional type will always go extinct, and in which the PFT destined for extinction is pre-determined by the parameterization and is independent of the initial conditions. In these cases, there is a fundamental dichotomy between what the model equations suggest is the case (in this case a system with two competing phytoplankton) and how the model actually functions (as a single phytoplankton–zooplankton predator–prey system). Two non-PFT systems, in which the real eigenvector of one of the predator–prey critical points on the faces is positive and the other one is negative, are shown in Figs 6 and 7. In these cases, there is no predator–prey–prey critical point in the feasible region of the state space and the phytoplankton functional type model is illusory, in that more than one PFT exists only in the initial transient dynamics; ultimately one PFT will always dominate and the other PFT will always go extinct. In contrast to the pseudo-PFT systems, extinction in these situations is determined not by the initial populations, but by an intrinsic property of the parameter set chosen, and there are no circumstances under which the loser can survive. In Fig. 6, for example, the predator–prey critical point on the (P₁, Z) face has a negative real eigenvalue while the predator–prey critical point on the (P₂, Z) face has a positive real eigenvalue. All initial conditions are therefore eventually attracted to the stable spiral critical point on the (P₁, Z) face ("extinction" of P₂), even though the initial dynamics of the rightmost trajectory suggest that P₂ may dominate.

Figure 7 shows the reverse situation, when the stable spiral on the (P₁, Z) face has a positive real eigenvalue and the stable limit cycle on the (P₂, Z) face has a negative real eigenvalue. In this case, all initial conditions result in the "extinction" of P₁. The dynamics in Figs 6 and 7 are consistent with the "Competitive Exclusion Principle" first articulated by Gause (Gause, 1934).

Sensitivity of the predator–prey–prey critical point

The results of the numerical investigation of the sensitivity of the predator–prey–prey critical point for the bona fide PFT system of Fig. 2 are summarized in Table V. The minimum and maximum values in Table V refer to the values of each parameter for which the predator–prey–prey critical point (E) leaves the ecologically feasible region of the state space. These are defined analytically by equation (A.29). The effect of parameter variations on the predator–prey–prey critical point is consistent across all parameters and in accord with the critical point and eigenvalue analysis above. The sensitivity analysis identifies the range of parameter values for which E exists in the interior of the parameter space.

View this table: [in this window] [in a new window]   Table IV: Criteria for classification of parameter sets for the NP1P2Z system

 

View this table: [in this window] [in a new window]   Table V: Summary of sensitivity analysis of the interior predator–prey–prey (E) critical point of the NP1P2Z model based on the parameter values used in Fig. 2 (column 1 Table II)

 
Variations in the parameter values move the predator–prey–prey critical point (E) across the state space between the predator–prey critical points (D and F). The parameter values for which E is located within the feasible region of the state space are listed in Table V, with the sensitivity metric where available and a ranking of relative sensitivities to one-at-a-time parameter variations.

We observe that because we have used a parameter set that results in a stable spiral on one face and a limit cycle on the other, the sensitivity of the critical point E is intimately related to the bifurcation behaviour of the system. As E moves across the state space, it undergoes a Hopf bifurcation near the centre of the parameter space. A Hopf bifurcation occurs where the real parts of the complex eigenvalues pass smoothly through zero. Figure 2 shows a system where E lies on the stable spiral side of the Hopf bifurcation while Fig. 3 shows a system where E lies on the unstable spiral side of the Hopf bifurcation.

We further note that when the critical point E leaves the ecologically feasible region of the state space in response to variations in the parameter values, it always leaves either via the predator–prey critical point (D) on the (P₁, Z) face or via the predator–prey critical point (F) on the (P₂, Z) face. Mathematically, this is referred to as a transcritical bifurcation. At a transcritical bifurcation, the interior predator–prey–prey critical point occupies the same position as the predator–prey critical point on the face, and the two critical points exchange stability (i.e. the real parts of their eigenvalues simultaneously change sign in opposite directions). Therefore, after the transcritical bifurcation, the predator–prey critical point on the face is stable in the (P₁, P₂) dimension and the dynamics of the system are confined to the face. The predator–prey–prey critical point (E) now lies outside the ecologically feasible region of the state space and has no influence on the dynamics of the system.

The sensitivity analysis of the bona fide PFT system further emphasizes its delicately balanced nature, as subtle variations in most parameters can move the predator–prey–prey critical point through the Hopf bifurcation and across the entire width of the state space in the (P₁, P₂) dimension. One-at-a-time variations in parameter values fail to do this only when a parameter loses ecological realism (becomes negative or, in the cases of and , larger than one) prior to the critical point reaching the face. We note that our analysis does not exclude the possibility that a robust parameter set, where parameters can be varied dramatically without affecting the dynamics, could exist for this model. However, our results indicate that the existence of such a parameter set is unlikely.

Ubiquity of the predator–prey–prey critical point

The criteria for determining valid parameter sets for the NP₁P₂Z system, derived from the analysis presented in the Appendix, are presented in Table IV. The criteria for determining those parameter sets that will result in bona fide PFT dynamics are listed in the lower part of Table IV. The 5 x 10⁶ parameter sets randomly generated in the search of the parameter space returned 1.9 x 10⁶ (about 38%) invalid parameter sets because they did not meet Kolmogorov's criteria of producing an autotroph and a predator–prey critical point in both of the NP₁Z and NP₂Z systems (Table IV). The distribution of the eigenvalues _{D –2} and _{F –1} for the 3.1 x 10⁶ (62%) of parameter sets that were valid is shown in Fig. 8.

View larger version (29K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 8. Two-dimensional frequency distribution of the values of the eigenvalues D –2 and F –1 randomly sampled from the parameter space defined for the sensitivity analysis (see text for details).

 
Eigenvalues associated with bona fide PFT dynamics lie in the positive–positive quadrant at the front of Fig. 8. These are the 17 767 (0.6%) valid parameter sets that had both the eigenvalues _{D –2} and _{F –1} positive, indicating that the system had a stable interior predator–prey–prey critical point and would exhibit bona fide PFT system behaviour such as in Figs 2 and 3. Pseudo-PFT systems were also relatively rare with only a further 228 389 (7.4%) parameter sets having both eigenvalues negative (rear quadrant of Fig. 8), indicating that these systems had a horizontally unstable interior predator–prey–prey critical point. These systems would exhibit the "pseudo"-PFT behaviour of Figs 5 and 6, where one or the other functional type will always go extinct depending on the initial conditions.

The majority of parameter sets lie in the positive–negative (left and right) quadrants of Fig. 8. These parameter sets, which had one positive and one negative eigenvalue, dominated the sample, representing 92% of valid parameter sets. These parameter sets do not result in realistic PFT systems, as when implemented with these values the NP₁P₂Z system does not have an interior predator–prey–prey critical point (E). In these cases, one P always goes extinct for all initial conditions, with the P destined for extinction being predetermined by the parameters. This analysis reveals that the majority of possible parameter sets are degenerate in the context of PFT modelling.

Figure 8 clearly reveals that the eigenvalues are not distributed uniformly, but are concentrated near zero and in rows where one eigenvalue is small and positive and the other is negative with a larger range of magnitudes. NP₁P₂Z systems that use parameter sets with _{D –2} and _{F –1} near zero will be sensitive in the context that they lack strong determinants of their dynamics. The dynamics of these systems may be easily influenced by external factors, and they may have transient dynamics that give the appearance of bona fide PFT system behaviour. Such systems may require long and accurate numerical integrations to reveal their true nature, but the outcomes of these integrations may also be critically affected by variations in the physical environment in which the model ecosystem is embedded.

The fact that if _{D –2} and _{F –1} have opposite signs, they will usually comprise a large magnitude negative eigenvalue and a small magnitude positive eigenvalue reveals that competition is finely balanced. Whether the eventual winner can grow or not is essentially determined by the nature of the loser's critical point. If the loser and its predators have characteristics that leave significant levels of inorganic nutrient available for a competitor, and does not support a substantial predator population that could also graze on the competitor, then the competitor has an opportunity to grow even from almost zero populations. In these cases, once the competitor capitalizes on its opportunity to grow, the competition is over and the whole system then rapidly moves to the winner's face along its strongly negative eigenvector.

The frequency distributions of the parameters that contributed bona fide PFT systems are shown in Fig. 9. This reveals that bona fide PFT parameter values were generally found throughout the parameter spaces; no parameter combination has a zero frequency. However, bona fide PFT dynamics are difficult to find if some parameters are approximately equal, for example, the distribution of zooplankton grazing coefficients that revealed few parameter sets for which . Similarly, the distribution of phytoplankton half-saturation constants revealed a paucity of parameter sets for which .

View larger version (36K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 9. Two-dimensional frequency distributions of the parameter values of those sets randomly sampled from the parameter space defined for the sensitivity analysis (see text for details) that had both D –2 and F –1 positive. The frequency distributions are calculated for pairs of equivalent parameters for P1 and P2.

 
The distributions therefore reveal that, somewhat counter-intuitively, bona fide PFT systems are most likely to have phytoplankton with similar maximum growth rates (). However, multiple phytoplankton with similar growth rates appear sustainable only if they differ in other properties. Competing phytoplankton are more likely to coexist if they have quite different grazing susceptibilities (or conversely, the zooplankton must have a preference with ), and different strategies for utilizing nutrients ().

Bona fide PFT systems with both phytoplankton having low specific mortality rates () are also rare, although the PFT parameter sets are uniformly distributed throughout the rest of the parameter space. A similar distribution is observed for the zooplankton assimilation efficiencies (), although to a lesser degree; in this case a substantial number of PFT sets had equally low values for both and , and few parameter sets had very large values () of one or both.

The opposite distribution is observed for and , where most PFT systems had low values for both and , and few had high values for both. This suggests that systems in which both predator–prey interactions have the form of Lotka–Volterra grazing (i.e. stable spiral dynamics) are more likely to exhibit bona fide PFT dynamics than systems that both have Michalis–Menten grazing interactions (with stable limit cycle dynamics). Systems that had one Lotka–Volterra and one Michalis–Menten sub-system had intermediate likelihood, as did sub-systems with grazing interactions intermediate between solely Lotka–Volterra and solely Michalis–Menten. However, we note that the form of the frequency distribution of and values is a slightly tilted plane, and that therefore while spirally stable systems may be more likely to have coexisting competing species, it would not be appropriate to utilize this as a primary heuristic for constructing PFT models.

A consistent feature of all the plots in Fig. 9 is that none of the frequency distributions is zero at any point. This indicates that, while some parameter combinations are more likely to result in bona fide PFT dynamics than others, no unequivocal rules that will guarantee bona fide PFT dynamics exist. Conversely, parameter sets that result in bona fide PFT dynamics appear to exist throughout the parameter space, although more dense in some areas than others.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
The examination of the dynamics of a simple PFT model with two competing phytoplankton functional types has provided some useful insights into, and heuristics for, the construction of more complex PFT models that can be used to simulate biogeochemical cycling in the oceans. Our analysis reveals that bona fide PFT systems are rare in the ecological parameter space, and that it requires a finely tuned and delicately balanced parameter set to maintain both phytoplankton functional types extant in the model simulations. However, our results also provide practical heuristics that can allow model developers to check that the systems they devise from biological criteria have dynamical characteristics that maintain all phytoplankton extant. These heuristics apply to Kolmogorov systems (systems that can be written in the form ) and that conserve mass (u_i = constant) and we reiterate that many ecological models currently in use comply with these criteria. Most exceptions will occur because models are designed to explicitly represent losses of mass from the mixed layer due to export to the deep ocean; however, we suggest that the dynamics of these models may be assessed for the special case where exports are exactly balanced by imports of nutrient, and hence mass is conserved.

The application of Kolmogorov's (Kolmogorov, 1936) criteria that each predator–prey sub-system must have both a feasible autotroph critical point and a feasible predator–prey critical point provides a useful heuristic to constrain a search of ecological parameter spaces that otherwise have quite nebulous bounds. Similarly, the eigenvalues of the inward-pointing eigenvectors of the predator–prey critical points on the faces of the state space can provide a useful goal function in a search for parameter sets that result in bona fide PFT systems. Once a model has been developed from biological principles, bona fide PFT dynamics can in many cases be assured by first constraining the parameter space and then evaluating the growth functions (f_i) of the other phytoplankton at each of the critical points on the boundaries (edges or faces) of the state space. This relatively simple calculation gives the eigenvalues of the eigenvectors that are orthogonal to each predator-prey face of the state space. If every critical point has at least one of these eigenvalues (f_i) greater than zero, then the system will exhibit bona fide PFT dynamics, that is, no phytoplankton will go extinct. This heuristic can potentially be applied to any-dimension Kolmogorov system, as the eigenvalues of the PFTs that equal zero at a critical point on a face are given by their f_i evaluated at the critical point, and their eigenvectors are orthogonal to the face. However, difficulties in evaluating the critical points on the faces may impose constraints on application to some higher dimension systems.

The fact that the eigenvalues of the predator–prey critical points on the faces of the state space provide useful indications of the existence and stability of the interior predator–prey–prey critical point is interesting in a more general context. The possibility for a competing phytoplankton to bloom and hence destabilize a predator–prey critical point on a face is determined by how much resource its competitor leaves available to it, and how large a population of grazers it supports. If too little resource is available, or too many grazers abound (i.e. f_i < 0), then the potential competitor will not be able to bloom and will become extinct. A clear implication of this research is that extant ecosystems are not random assemblages of species and that consequently estimates of the parameter values that describe PFT traits need to be refined as systems.

The result that parameter sets that allow for continued co-existence of multiple phytoplankton functional types in PFT models are rare should perhaps not come as a surprise, as evolution has had many, many generations of plankton to operate on to achieve the optimization of the exploitation of niches that this implies. While it may be reasonable to suppose that the delicate balance of biological and ecological properties required to allow the coexistence of many species can be achieved by evolution, it is another matter for modellers to derive analogous parameter sets for their models. This has been the subject of recent discussion, with concerns being expressed that the demands of climate modelling for the implementation of PFT models has progressed in advance of their validation (Anderson, 2005; Anderson, 2006; Flynn, 2006; Le Quéré, 2006). Our results suggest that all the approaches suggested by these authors are relevant to the development of robust PFT models. Le Quéré (Le Quéré, 2006) observed parameter estimates for PFT models can be refined by laboratory and field measurement and used as the basis for model calibration. We concur with Flynn (Flynn, 2006) who considered that these measurements must be taken for groups of organisms. Our results reveal that the parameter spaces so defined can be expected to yield a number of finely tuned possible parameterizations that will result in a PFT model with bona fide dynamics; a first step in determining which of these parameter sets is best for simulations should be to rigorously validate the model outputs against observed data, as suggested by Anderson (Anderson, 2006).

Our work also provides some insights into how PFT models could be made more robustly persistent, such that parameter sets that result in bona fide PFT dynamics are easier to find. One obvious option is the use of nonlinear mortality terms. The linear mortality terms used in the NP₁P₂Z model result in a negative constant in each f_i, whereas the use of nonlinear mortality terms would result in the mortality constants appearing in the f_i as multipliers of the state variables. As we are interested in the f_i of the state variables where they equal zero, these terms vanish at the relevant critical points and allow greater freedom in the values of the growth and grazing parameters. Clearly, given the large literature on the effect of quadratic zooplankton mortality terms on the dynamics of ecological models (Steele and Henderson, 1992; Edwards and Brindley, 1999; Edwards and Yool, 2000), the use of nonlinear phytoplankton mortality terms may represent a non-trivial change to the dynamics of the system. However, we observe that quadratic mortality terms have allowed better representation of phytoplankton dynamics in complex models (Raick et al., 2006), and have been demonstrated to prevent competitive exclusion in simple models (Ruan et al., 2007). Our preliminary investigations revealed that the inclusion of quadratic phytoplankton mortality terms in this model more than doubled the number of bona fide parameter sets.

Hutchinson (Hutchinson, 1961) considered that environmental fluctuations may allow coexistence of species that would not occur in stable environments, and this has been demonstrated in some models (Scheffer et al., 2003). The principal environmental forcing in plankton ecosystems is the seasonal change in the irradiance field experienced by phytoplankton due to changes in solar irradiance at the sea surface and changes in the depth of the mixed layer. This serves to change the maximum phytoplankton growth rate, and may affect the signs of the eigenvalues of competing phytoplankton (i.e. the sign of the f_i on the opposite predator–prey face). Similarly, changes in parameter values due to variations in assimilation efficiency and food quality that are known to occur in plankton ecosystems (Mitra et al., 2007) may affect the persistence of some PFTs as these also affect the values of the f_i that determine persistence or extinction.

Mathematical analysis of the effects of seasonal changes in the eigenvalues associated with the competing phytoplankton is not straightforward. However, we could expect that such changes could significantly affect the co-existence properties of PFT systems if the eigenvalues of competing phytoplankton changed sign in a regular way. If a small positive eigenvalue became predominately negative as a result of seasonal forcing, then the associated PFT may become extinct, and vice versa. Seasonal forcing could therefore have potentially the most influence on pseudo-PFT systems, as analysis of the NP₁P₂Z model suggests that non-PFT systems have eigenvalues that typically occur in pairs with one positive and small and the other negative and large. It is less likely that a seasonal forcing would change a strongly negative eigenvalue than a weakly negative one. Hence pseudo-PFT systems, which analysis of the NP₁P₂Z model suggests typically have small, negative eigenvalues, may be most likely to have qualitative dynamics susceptible to external forcing. Seasonal forcing of these systems may alter the eigenvalues of the predator–prey points to become mostly positive, and preserve both phytoplankton. Therefore, some systems that do not exhibit endogenous bona fide PFT dynamics, when subjected to particular physical forcings, may retain phytoplankton that would go extinct in unforced scenarios.

The availability of sophisticated nonlinear optimization tools allows modellers to find by brute force parameter sets that result in bona fide PFT models that reproduce observed dynamics. However, it remains a moot point whether PFT models parameterized in such way can respond appropriately in scenarios of a changing environment, such as proposed under global climate change. In this context, an appropriate response means that the eigenvalues of the interior critical point and the corresponding critical points on the faces change sign appropriately in response to changes in the environment to reflect extinctions and changes in the community structure.

Finally, we reiterate that we have only considered the case of two competing autotrophs with a common grazer in our analysis of the generic three-dimensional Kolmogorov system. Other generic three-dimensional Kolmogorov systems that warrant investigation are systems with one autotroph that has two competing grazers (i.e. an NPZ₁Z₂ system), and systems with one grazer and an autotroph that is potentially limited by one of two nutrients (an N₁N₂PZ system). These cases will be considered in forthcoming work. There are, of course, many higher dimension PFT models in existence, and a theme of future work will be to investigate the extent to which the techniques we have used in this work can be extended to higher dimensions.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
The dynamical systems perspective that we have taken in this paper has provided insights into how to construct more robust PFT models for application in modelling climatically important biogeochemical cycles in the ocean. Our analysis of the generic Kolmogorov system has generated some useful generic heuristics to guide modellers when developing and parameterizing PFT models. These heuristics include:

• Models that explicitly resolve nutrients, as all biogeochemical models do, must conserve nutrients in order to be written in Kolmogorov form. In many cases this is a reasonable assumption, as plankton ecosystems recycle nutrients very tightly within the mixed layer. However, in models that seek to explicitly resolve movements of nutrients between the surface and the deep ocean, it is suggested that the special case where imports and exports are balanced be analysed.
  
• Complicated PFT models may be broken down into simpler predator–prey and competition sub-models. These sub-systems control the dynamics on the vertices, edges and faces of the full system. Each sub-system of a Kolmogorov system is also a Kolmogorov system.
  
• The essence of Kolmogorov's criteria for reasonable predator–prey dynamics, that each system has critical points at the origin (everything extinct), an autotroph-only point, and a predator–prey point within the feasible state space, may be used to constrain the parameter space to be searched.
  
• The eigenvalues of the "competition" eigenvectors that are orthogonal to each predator–prey face provide a useful indicator of the existence and stability properties of the interior point that is fundamental to a bona fide PFT model. Analytic expressions for these eigenvalues are very easily found in Kolmogorov systems, being given by the f_i of the competitor evaluated at the predator–prey critical point.
  
• Using nonlinear mortality terms for PFTs involved in competition may assist in making PFT systems more robust, in that nonlinear mortality may increase the number of potential parameter sets that result in bona fide PFT dynamics.
  

The application of our analysis to an example N₁N₂PZ Kolmogorov system revealed further heuristics that may also be useful for PFT modelling:

• Parameter sets that result in bona fide PFT systems with competing phytoplankton are rare: we found an order of magnitude more parameter sets that resulted in pseudo-PFT systems, and an order of magnitude again more parameter sets that resulted in non-PFT systems
  
• Bona fide PFT parameter sets are distributed throughout the parameter space; therefore, for any PFT model there may be many bona fide parameter sets with significantly different properties
  
• PFT models require finely tuned parameter sets in order to have bona fide PFT properties; a consequence of this is that near any bona fide parameter set there are many parameter sets that do not result in bona fide PFT dynamics
  

The ubiquity of parameter sets throughout parameter space presents a substantial challenge to PFT modelling in determining how PFT models should be parameterised. Choosing parameter sets from within measured ranges so that PFT models fit contemporary data, in a calibration and validation sense, is an obvious starting point. However, PFT models that are used in simulations of climate change, one of the primary drivers of the development of PFT models, must also respond appropriately to changes in their external environment. Our results suggest that the development of PFT models that can simulate the changes in community composition that this implies need, in addition to refined measurements of parameters and acquisition of longer time series of data, intimate knowledge of the inherent properties of the models.

We wish, in conclusion, to draw attention to the import of this work in the context of the current practice in biogeochemical modelling of plankton in the oceans. In the words of an anonymous reviewer of this work:

"... current PFT models have concentration thresholds below which PFTs are forced to survive by setting their loss terms to zero or even inverting the signs of the losses."

We would argue that such artifices may compromise the capacity of PFT models to simulate the feedback effect that plankton ecosystems may have on climate in response to climate change. These artifices do not validly represent real processes such as encystment, as is sometimes argued, because recovery from encystment still requires favourable growth conditions (f_i > 0) that, if they occurred, would render the artifice unnecessary. The recent results of Montes-Hugo et al. (Montes-Hugo et al., 2009) suggest that changes in plankton community structure are already occurring in response to climate change. It is unlikely that a PFT model that prevents PFTs from going extinct by enforcing minimum concentrations in the numerical solution may either simulate or predict such changes. The work we have presented here provides some indications of how we might build PFT models that obviate the requirement for such measures.

	SUPPLEMENTARY DATA

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
Supplementary data can be found online at http://plankt.oxfordjournals.org.

	FUNDING

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
The authors wish to thank the Zilkha Trust, Lincoln College, Oxford for supporting this research.

	APPENDIX

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
This material provides analytic solutions, as far as possible, for the critical points and their associated eigenvalues and eigenvectors for the example NP₁P₂Z system (equations (4)–(7)). Each critical point is labelled according to Fig. 1, and the eigenvalues at each critical point are additionally differentiated by numbering from 1 to 3.

ORIGIN CRITICAL POINT (O)

The eigenvalues of the NP₁P₂Z system at the origin are:

(A.1)

(A.2)

(A.3)

PREY-ONLY CRITICAL POINTS (A, C)

The NP₁P₂Z system has an autotroph critical point in each of the NP₁ (A) and NP₂ (C) sub-systems, on the P₁ and P₂ axes, respectively, at

(A.4)

and

(A.5)

Note that for parameter values such that or , or or , these points will lie outside the ecologically feasible region and in this case the systems no longer exhibit realistic dynamics. The eigenvalues of these points are given by, for the NP₁ sub-system (A):

(A.6)

(A.7)

(A.8)

and for the NP₂ sub-system (C):

(A.9)

(A.10)

(A.11)

DUAL PREY-ONLY CRITICAL POINT (B)

The NP₁P₂Z system has a dual autotroph critical point in the NP₁P₂ (B) sub-system, in the (P₁, P₂) plane, at

(A.12)

We note that consistency requires that In this case an autotroph critical line, every point of which satisfies equation (A.12), exists across the (P₁, P₂) dimension and intersects the autotroph critical points on each of the P₁ and P₂ axes:

(A.13)

This critical line has the eigenvalues

(A.14)

(A.15)

(A.16)

This reveals that the eigenvalue leading from the origin to the dual-autotroph critical point (_{B– 1}) is effectively a weighted average of the autotroph eigenvalues of the individual predator–prey critical points on the faces, and similarly the eigenvalue associated with the predator (_{B– 3}) is a weighted average of the predator eigenvalues of the individual critical points. The eigenvalue pointing across the (P₁, P₂) dimension between the two-autotroph points (_{B– 2}) is identically zero, indicating neutral stability in this direction. This is consistent with equation (A.13), which reveals that the line conserves mass with no distinction between . We note that when parameter sets have the final relative proportions of will be determined solely by their relative starting proportions (i.e. their initial conditions).

PREDATOR–PREY CRITICAL POINTS (D, F)

The NP₁P₂Z system has a predator–prey critical point on each of the (P₁, Z) and (P₂, Z) faces. The predator–prey critical point NP₁Z on the (P₁, Z) face (D) is located at

(A.17)

(A.18)

(A.19)

(A.20)

The NP₂ Z point on the (P₂, Z) face (F) is located at

(A.21)

(A.22)

(A.23)

(A.24)

The eigenvalues of these points are

(A.25)

(A.26)

and

(A.27)

(A.28)

respectively. The complex eigenvalues control the dynamics on the (P₁, Z) and (P₂, Z) planes, and are composed only of parameters relating to the autotroph and Z (i.e. the complex eigenvalues of the NP₁Z point include parameters that appear only in the P₁ and Z equations and vice versa).

PREDATOR–PREY–PREY CRITICAL POINT (E)

The NP₁P₂Z system has a predator–prey–prey critical point (E) that always lies on a line passing through the predator–prey critical points on the faces (NP₁Z and NP₂Z). This point may lie inside or outside the feasible region, depending on the parameter values chosen. The location of this predator–prey–prey critical point may only be obtained numerically by iteratively solving the following equations (A.29)–(A.31):

(A.29)

(A.30)

where

(A.31)

where

Note that equation (A.29) describes a relation between , and that when and when These conditions mark the endpoints of the range of values for which E is ecologically realistic.

	ACKNOWLEDGEMENTS

 
The authors would like to thank three anonymous reviewers whose commentary and advice led to substantial improvement in this manuscript.

	Notes

 
Corresponding editor: Roger Harris

	REFERENCES

TOP ABSTRACT INTRODUCTION METHOD RESULTS DISCUSSION CONCLUSIONS SUPPLEMENTARY DATA FUNDING APPENDIX REFERENCES

 
Anderson T. R. Plankton functional type modelling: running before we can walk? J. Plankton Res. (2005) 27:1073–1081.[Abstract/Free Full Text]

Anderson T. R. Confronting complexity: reply to Le Quéré and Flynn. J. Plankton Res. (2006) 28:877–878.[Free Full Text]

Billen G., Becquevort S. Phytoplankton–bacteria relationships in the Antarctic marine ecosystem. Polar Res. (1991) 10:245–253.[CrossRef][Web of Science]

Cropp R. A., Norbury J., Gabric A. J., et al. Modeling dimethylsulphide production in the upper ocean. Global Biogeochem. Cycles (2004) 18:1–21.

Dhooge A., Govaerts W., Kuznetsov Y. A. MATCONT: a Matlab package for numerical bifurcation analysis of ODEs. ACM Trans. Math. Soft. (2003) 29:141–164.[CrossRef]

Edwards A. M., Brindley J. Zooplankton mortality and the dynamical behaviour of plankton population models. Bull. Math. Biol. (1999) 61:303–339.[CrossRef][Web of Science][Medline]

Edwards A. M., Yool A. The role of higher predation in plankton population models. J. Plankton Res. (2000) 22:1085–1112.[Abstract/Free Full Text]

Falkowski P., Laws E. A., Barber R. T., et al. Phytoplankton and their role in primary, new and export production. In: Ocean Biogeochemistry.—Fasham M. J. R., ed. (2003) New York: Springer. 99–121.

Fenchel T. Ecology of heterotrophic microflagellates. IV. Quantitative occurrence and importance as bacterial consumers. Mar. Ecol. Progress Ser. (1982) 9:35–42.[CrossRef]

Flynn K. J. Modelling multi-nutrient interactions in phytoplankton: balancing simplicity and realism. Prog. Oceanogr. (2003) 56:249–279.[CrossRef]

Flynn K. J. Reply to Horizons article ‘Plankton functional type modelling: running before we can walk’ Anderson (2005): II. Putting trophic functionality into plankton functional types. J. Plankton Res. (2006) 28:873–875.[Free Full Text]

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

Gabric A. J., Matrai P. A., Vernet M. Modelling the production and cycling of dimethylsulphide during the vernal bloom in the Barents Sea. Tellus (1999) 51B:919–937.

Gabric A. J., Whetton P. H., Cropp R. A. Dimethylsulphide production in the subantarctic southern ocean under enhanced greenhouse conditions. Tellus (2001) 53:273–287.

Gause G. F. The Struggle for Existence (1934) Baltimore: Williams and Wilkins.

Hansen B., Christiansen S., Pedersen G. Plankton dynamics in the marginal ice zone of the central Barents Sea during spring: carbon flow and structure of the grazer food chain. Polar Biol. (1996) 16:115–128.[CrossRef][Web of Science]

Hood R. R., Laws E. A., Armstrong R. A., et al. Pelagic functional group modeling: progress, challenges and prospects. Deep-Sea Res. II (2006) 53:459–512.[CrossRef]

Huang X. C., Zhu L. Limit cycles in a general Kolmogorov model. Nonlinear Anal. (2005) 60:1393–1414.[CrossRef]

Hutchinson G. E. The paradox of the plankton. Am. Nat. (1961) 95:137–145.[CrossRef]

Kolmogorov A. N. Sulla Teoria di Volterra della Lotta per l'Esisttenza. Giorn. Instituto Ital. Attuari (1936) 7:74–80.

Le Quéré C. Reply to horizons article ‘phytoplankton functional type modelling: running before we can walk’ Anderson 2005: I – abrupt changes in marine ecosystems? J. Plankton Res. (2006) 28:871–872. doi:810.1093/plankt/fbl1014.[Free Full Text]

Le Quéré C., Harrison S. P., Prentice I. C., et al. Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models. Global Change Biol. (2005) 11:2016–2040.

May R. M. Stability and Complexity in Model Ecosystems (1973) Princeton: Princeton University Press.

Mitra A. Are closure terms appropriate or necessary descriptors of zooplankton loss in nutrient-phytoplankton-zooplankton type models? Ecol. Model. (2009) 220:611–620.[CrossRef][Web of Science]

Mitra A., Flynn K. J., Fasham M. J. R. Accounting for grazing dynamics in nitrogen–phytoplankton–zooplankton models. Limnol. Oceanogr. (2007) 52:649–661.

Moloney C. L., Bergh M. O., Field J. G., et al. The effect of sedimentation and microbial nitrogen regeneration in a plankton community: a simulation investigation. J. Plankton Res. (1986) 8:427–445.[Abstract/Free Full Text]

Montes-Hugo M., Doney S., Ducklow H. W., et al. Recent changes in phytoplankton communities associated with rapid regional climate change along the Western Antarctic Peninsula. Science (2009) 323:1470–1473.[Abstract/Free Full Text]

Muller-Niklas G., Herndl G. J. Dynamics of bacterioplankton during a phytoplankton bloom in the high Arctic waters of the Franz-Joseph Land archipelago. Aquat. Microb. Ecol. (1996) 11:111–118.[CrossRef]

Raick C., Soetaert K., Gregoire M. Model complexity and performance: how far can we simplify? Prog. Oceanogr. (2006) 70:27–57.[CrossRef]

Ruan S., Ardito A., Ricciardi P., et al. Coexistence in competition models with density-dependent mortality. C. R. Biol. (2007) 330:845–854.[CrossRef][Web of Science][Medline]

Scheffer M., Rinaldi S., Huisman J., et al. Why plankton communities have no equilibrium: solutions to the paradox. Hydrobiologia (2003) 491:9–18.[CrossRef][Web of Science]

Slagstad D., Stole-Hansen K. Dynamics of plankton growth in the Barents Sea. Polar Res. (1991) 10:173–186.[CrossRef][Web of Science]

Spitz Y. H., Moisan J. R., Abbott M. R. Configuring an ecosystem model using data from the Bermuda Atlantic Time Series (BATS). Deep-Sea Res. II (2001) 48:1733–1768.[CrossRef]

Steele J. H., Henderson E. W. The role of predation in plankton models. J. Plankton Res. (1992) 14:157–172.[Abstract/Free Full Text]

Vallina S. M., Simo R., Popova E. E., et al. A dynamic model of ocean sulfur (DMOS) applied to the Sargasso Sea: simulating the dimethylsulfide (DMS) summer paradox. J. Geophys. Res. (Biogeosciences) (2008) 113:1–23.

CiteULike    Connotea    Del.icio.us    What's this?

This Article Abstract FREE Full Text (PDF) Supplementary Data All Versions of this Article: 31/9/939    most recent fbp042v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Google Scholar Articles by Cropp, R. Articles by Norbury, J. PubMed Articles by Cropp, R. Articles by Norbury, J. Social Bookmarking          What's this?

Online ISSN 1464-3774 - Print ISSN 0142-7873

Copyright © 2009 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals China Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics