INTRODUCTIONTop
A feature of phytoplankton communities is the occurrence of population explosions or blooms. Such blooms may be an integral component of marine ecosystem dynamics, but sometimes they can have negative impacts on other marine organisms, on human health or on the socioeconomic wellbeing of the area in which the bloom occurs. Basic physicochemical conditions for bloom development include nutrient availability and low physical dispersion losses (Kierstead and Slobodkin 1953Kierstead H., Slobodkin L. 1953. The size of water masses containing plankton blooms. J. Mar. Res. 12: 141147., Wyatt and Horwood 1973Wyatt T., Horwood J. 1973. Model which generates red tides. Nature 244: 238240., Margalef et al. 1979Margalef R., Estrada M., Blasco D. 1979. Functional morphology of organisms involved in red tides, as adapted do cecaying turbulence. In: Taylor D.L. and Sliger H.H. (eds), Toxic dinoagellate blooms. Elsevier, North Holland, pp. 315320., Smayda and Reynolds 2001Smayda T.J., Reynolds C.S. 2001. Community assembly in marine phytoplankton: application of recent models to harmful dinoagellate blooms. J. Plankton Res. 23(5): 447461.). These factors interact with community interactions to control the net population growth of a species. In particular, differential grazing of phytoplankton species by herbivores may play a key role in regulating the abundance and composition patterns in the phytoplankton community (Steele 1974Steele J. 1974. The structure of marine ecosystems. Harvard University Press., Holt et al. 1994Holt R., Grover J., Tilman D. 1994. Simple rules for interspecific dominance in systems with exploitative and apparent competition. Am. Nat. 144(5): 741771., Leibold 1996Leibold M.A. 1996. A graphical model of keystone predators in food webs: trophic regulation of abundance, incidence, and diversity patterns in communities. Am. Nat. 147: 784812., Armstrong 2003Armstrong R.A. 2003. A hybrid spectral representation of phytoplankton growth and zooplancton response: The “control rod” model of plankton interaction. Deep Sea Res. Part II 50: 28952916.). Prey selection may allow or inhibit the proliferation of certain phytoplankton species, as has been shown for high nutrientlow chlorophyll (HNLC) regions (Leising et al. 2005Leising A., Horner R., Pierson J., et al. 2005. The balance between microzooplankton grazing and phytoplankton growth in a highly productive estuarine fjord. Progr. Oceanogr. 67: 366383.). Some harmful algal blooms (HABs) have been attributed to grazing impairment due to toxin production. Wellknown examples include the Chrysochromulina polylepis red tide studied by Maestrini and Granéli (1991)Maestrini S., Granéli E. 1991. Environmental conditions and ecophysiological mechanisms which led to the 1988 Chrysochromulina polylepis bloom: an hypothesis. Ocean. Acta 14: 397413. and the persistent brown tides in Texas of the chrysophyte microalga Aureoumbra lagunensis, which appears to be toxic for some grazers and of poor quality as food for others (Buskey and Hyatt 1995Buskey E., Hyatt C. 1995. Effects of the Texas (USA) ‘brown tide’ alga on planktonic grazers. Mar. Ecol. Prog. Ser. 126: 285292.). Some blooming species may escape control by microzooplankton through a combination of grazing avoidance mechanisms (such as larger size, colonies, spines and toxic compounds) and factors such as poor food quality (Irigoien et al. 2005Irigoien X., Flynn K.J., Harris R.P. 2005. Phytoplankton blooms: a loophole in microzooplankton grazing impact? J. Plankton Res. 27: 313321.).
An interesting situation arises when one among two otherwise similar microalgal competitors has a lower intrinsic growth rate than the other, but possesses a defence mechanism that affects prey selectivity by grazers or predators, inducing a decreased grazing pressure for itself and an increased grazing pressure for the competitor. Such tradeoffs can have important impacts on food web dynamics and have been the subject of numerous studies (Leibold 1989Leibold M.A. 1989. Resource edibility and the effects of predators and productivity onTillmann U., Hesse K., Colijn F. 2000. Planktonic primary production in the German Wadden Sea. J. Plankton Res. 22: 12531276., Jessup and Bohannan 2008Jessup C.M., Bohannan B.J.M. 2008. The shape of an ecological tradeoff varies with environment. Ecol. Let. 11: 947959.). Tradeoffs can arise, for example, if toxin production reduces the amount of energy or nutrients allocated to growth, or if a lower nutritional quality of one of the competitors is linked with reduced growth rate (Yoshida et al. 2004Yoshida T., Hairston N.H., Ellner S. 2004. Evolutionary tradeoff between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris. Proc. R. Soc. London 271: 19471953., Litchman and Klausmeier 2008Litchman E., Klausmeier C. 2008. Traitbased community ecology of phytoplankton. Annu. Rev. Ecol. Evol. Syst. 39: 615639.). Another example would be the metabolic costs associated with the production of morphological defences (van Donk 1997van Donk E. 1997. Defenses in phytoplankton against grazing induced by nutrient limitation, UVB stress and infochemicals. Aq. Ecol. 31: 5358.). The importance of ’prey stoichiometry’ in predatorprey interactions has been well documented (e.g. Flynn et al. 1996Flynn K.J., Davidson K., Cunningham A. 1996. Prey selection and rejection by a microflagellate; implications for the study and operation of microbial food webs. J. Exp. Mar. Biol. Ecol. 196: 357372., Mitra and Flynn 2005Mitra A., Flynn K. 2005. Predatorprey interactions: is ‘ecological stoichiometry’ sufficient when good food goes bad? J. Plankton Res. 27: 393399., 2006Mitra A., Flynn K. 2006. Promotion of harmful algal blooms by zooplankton predatory activity. Biol. Lett. 2: 194197., Flynn 2010Flynn K.J. 2010. Do external resource ratios matter? Implications for modelling eutrophication events and controlling harmful algal blooms. J. Mar. Syst. 83: 170180.). Nutrient stress may also interact with growth rate and toxin production (Flynn 2002Flynn K.J. 2002. Toxin production in migrating dinoagellates: a modelling study of PSP producing alexandrium. Harmful Algae, 1: 147155., Flynn and Iringoien 2009Flynn K.J., Iringoien X. 2009. Why aldehydeinduced insidious effects cannot be considered as a diatom defence mechanism against copepods. Mar. Ecol. Prog. Ser. 377: 7989., Touzet et al. 2007Touzet N., Franco J., Raine R. 2007. Influence of inorganic nutrition on growth and PSP toxin production of Alexandrium minutum (Dinophyceae) from Cork Harbour, Ireland. Toxicon 50: 106119.). The role of grazingdeterrent toxins in allowing a toxic algal species to develop a bloom in competition with a nondefended competitor with higher growth rate was explored in a general framework by Solé et al. (2006b)Solé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664., by means of a simple predatorprey model. These authors used a logistic function for algal growth and a Holling II–type grazing response (Holling 1959Holling C. 1959. Some characteristics of simple types of predation and parasitism. Can. Entom. 91: 385398.) and found that the dominance of the toxic species could be predicted by a theoretical criterion that included a simple prey availability parameter indicating the degree of feeding selectivity by grazers.
In this work, we extend this approach to other functional response expressions, with the aim of establishing the parameter ranges that would allow an algal species with certain grazing defence mechanisms to dominate a nondefended competitor with a higher intrinsic growth rate. This assessment of the implications of functional response expressions and parameter choices may be helpful in the implementation of mathematical models of microalgal blooms.
Our approach is not dependent on the specific type of defence, which, as mentioned above, can be expressed in many ways, including reduced palatability. As discussed by Gentleman et al. (2003)Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875., Sailley et al. (2013)Sailley S., Vogt M., Doney S., et al. 2013. Comparing food web structures and dynamics across a suite of global marine ecosystem models. Ecol. Modell. 261262: 4357. and Cropp and Norbury (2013)Cropp R., Norbury J. 2013. Modelling plankton ecosystems and the Library of Lotka. In: Blackford J., Allen I., et al., Advances in Marine Ecosystem Modelling Research (III). J. Mar. Syst. 125: 313., ecosystem models may be very sensitive to the particular formulations adopted to express the functional response of zooplankton, which often imply untested assumptions. For example, Leising et al. (2005)Leising A., Horner R., Pierson J., et al. 2005. The balance between microzooplankton grazing and phytoplankton growth in a highly productive estuarine fjord. Progr. Oceanogr. 67: 366383. reported constraints of the microzooplankton grazing response functions that were necessary to achieve a HNLC situation. We will first consider a simplified system, composed of two phytoplankton prey and one zooplankton grazer species in a prebloom situation using general grazing and growth functions. Then we will explore the constraints that these functions must satisfy in order to allow one of the phytoplankton species to eventually become dominant. For simplicity, nutrient limitation will be considered as implicit in the value chosen for phytoplankton growth rates (Solé et al. 2006aSolé J., Estrada M., GarcíaLadona E. 2006a. Biological controls of Harmful Algal Blooms: A modelling study. J. Mar. Syst. 61: 165179.,bSolé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664.). We will examine the case of a species with a lower intrinsic growth rate than a competitor but with certain defence mechanisms that can reduce grazing pressure on it. In the following section we will apply these constraints to selected grazing functions, numerically solving the simple threespecies model to illustrate how changes in grazing control may trigger a bloom. To explore the applicability of the deduced relationships in an ecosystem context, we will carry out numerical simulations with the ERSEM model. Finally, we will discuss the obtained results and their implications for the modelling of algal blooms.
MATERIALS AND METHODSTop
Because our main interest is prebloom conditions, we will focus on the short time response starting from a nearequilibrium point (sources balanced by losses) and will examine the parameter constraints imposed on the grazing function if a particular phytoplankton species (the ‘dominant’ species) must attain a higher net growth rate than another (the ‘nondominant’ species) with a higher intrinsic growth rate (of course, the case is trivial if the ‘dominant’ species has a higher growth rate than the ‘nondominant’ one).
General conditions on feeding response functions
We propose a simplified ecosystem model that includes two microalgal populations, a nondominant (P_{1}) and a dominant one (P_{2}), and a zooplankton herbivore (Z) able to graze on both algal species. The general dynamics of this system can be formulated as in Solé et al. (2006b)Solé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664.:
$$\begin{array}{c}\frac{d{P}_{1}}{dt}={B}_{1}({P}_{1}){D}_{1}({P}_{1},{P}_{2},Z)\\ \frac{d{P}_{2}}{dt}={B}_{2}({P}_{2}){D}_{2}({P}_{1},{P}_{2},Z)\\ \frac{dZ}{dt}={D}_{1}({P}_{1},{P}_{2},Z)+{D}_{2}({P}_{1},{P}_{2},Z)E(Z)\end{array}$$ 
(1) 
B_{i} and D_{i} are, respectively, the functions describing the growth of the microalgae and the functional response of the zooplankton grazer. E(Z) is a function expressing zooplankton respiration and mortality. The B_{i} satisfy B_{i} (P_{i}=0)=0 and are increasing functions for small P_{i}. Similarly, the D_{i} satisfy D_{i} (P_{i} , P_{j} , Z=0)=0 and D_{i} (P_{i}=0, P_{j} , Z )=0 (grazing is null if there is either no prey or no grazer) and are also increasing functions of P_{i} and Z at least over a short period of time. For simplicity, in the system (1) any sloppy feeding or metabolic losses by zooplankton, as well as mortality due to predators or cannibalism, are assumed to be encapsulated in the E(Z) function. No competition terms are included because we assume that at this prebloom stage the nutrients do not have a limiting effect on algal growth.
In principle, all functions and parameters should be timedependent following physiological or environmental changes. However, as our study will address the shortterm response of the system, starting from prebloom conditions with low biomass concentration, we will assume that the system parameters remain constant at these prebloom conditions, before interspecific competition or environmental limitations become significant.
We consider that the prebloom system is near an equilibrium point (P_{1}^{*}, P_{2}^{*}, Z^{ }^{*}), in which source and loss rates are balanced so populations do not vary over time:
$$\left(\frac{d{P}_{1}}{dt}=\frac{d{P}_{2}}{dt}=\frac{d{Z}_{1}}{dt}=0\right)$$ 
(2) 
This state satisfies
B_{1} (P_{1}^{*}) = D_{1} (P_{1}^{*}, P_{2}^{*}, Z ^{*}) 
(2a) 
B_{2} (P_{2}^{*}) = D2 (P_{1}^{*}, P_{2}^{*}, Z ^{*}) 
(2b) 
E(Z^{*}) = D_{1} (P_{1}^{*}, P_{2}^{*}, Z ^{*}) + D_{2}(P_{1}^{*}, P_{2}^{*}, Z ^{*}) 
(2c) 
The following relationships can be derived from this system of equations near an equilibrium point:
 Relationship I: Taking the partial derivative with respect to P_{1} in (2c), we obtain
$$\frac{\partial {D}_{1}}{\partial {P}_{1}}=\frac{\partial {D}_{2}}{\partial {P}_{1}}$$ 
(3) 
which means that either
$$\frac{\partial {D}_{2}}{\partial {P}_{1}}=\frac{\partial {D}_{1}}{\partial {P}_{1}}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{or}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{1}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial {P}_{1}}$$ 
have an opposite sign. If
$$\frac{\partial {D}_{2}}{\partial {P}_{1}}=\frac{\partial {D}_{1}}{\partial {P}_{1}}=0$$, 
D_{1} and D_{2} do not depend on P_{1}, implying that grazers control only the growth of P_{2}, leaving P_{1} unbounded by grazing pressure. However we have assumed that both grazing functions depend on P_{1} and P_{1}; therefore,
$$\frac{\partial {D}_{1}}{\partial {P}_{1}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{1}}$$ 
must have the opposite sign and the second option holds.
Similarly, taking the partial derivative respect to P_{2}, we obtain
$$\frac{\partial {D}_{2}}{\partial {P}_{2}}=\frac{\partial {D}_{1}}{\partial {P}_{2}}$$ 
(4) 
Using the same argument as before, we deduce that
$$\frac{\partial {D}_{2}}{\partial {P}_{2}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial {P}_{2}}$$ 
must have the opposite sign. These results are in agreement with the properties of almost all grazing functions of Classes 1 and 2 studied by Gentleman et al. (2003)Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875..
– Relationship II: I holds for a system near an equilibrium point. However, we are also interested in how the system responds after a small perturbation. In this case, we can assume that the derivatives satisfy
$$\frac{{P}_{1}}{dt}\simeq 0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{{P}_{2}}{dt}\simeq 0$$ 
A necessary condition for the dominant species, in our case P_{2}, to bloom is to have a greater net growth rate than P_{1}
$$\frac{{P}_{1}}{dt}<\frac{{P}_{2}}{dt}\simeq 0$$ 
According to the system (1), we have
B_{1} − D_{1} < B_{2} − D_{2}, 
(5) 
Taking the derivatives (on P_{1} and P_{2}) of this expression and using I we can obtain the relationships:
$$\frac{d{B}_{2}}{d{P}_{2}}\text{\hspace{0.17em}}\text{>2}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{2}}$$ 
(6) 
$$\frac{d{B}_{1}}{d{P}_{1}}\text{\hspace{0.17em}}<\text{2}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial {P}_{1}}$$ 
(7) 
which bound the derivatives of the growth functions of the two phytoplankton species. Thus, determination of the growth functions and parameters implies a constraint on how the grazing functions must vary in order for one prey to dominate over the other after a small perturbation around the equilibrium state.
– Relationship III: From expressions (6) and (7), using Relationship I, we obtain:
$$\frac{\partial {D}_{1}}{\partial {P}_{1}}\text{\hspace{0.17em}}\text{>}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{2}}$$ 
(8) 
which implies that, near the equilibrium point, when a small perturbation is applied to the system, the grazing function D_{1} must grow faster in the P_{1} direction (the P_{1} axis in the phase space) than the grazing function D_{2} in the P_{2} direction.
These relationships are general and independent of the particular functions and parameters adopted. In the next section, we will examine whether the system is stable (it remains near equilibrium after a small perturbation) or unstable (it will diverge with time after a perturbation). This aspect is important in our study because blooms are effectively the result of a perturbation that separates the system from an equilibrium state and we want to know whether the populations (which we can assume to be initially near an equilibrium point) will evolve towards much higher values (a bloom situation) or remain bounded around zero.
Stability
The stability of the equilibrium solutions in (1) can be examined using phase plane analysis (Murray 1989Murray J. 1989. Mathematical Biology. Springer, Berlin., Truscott and Brindley 1994Truscott J., Brindley J. 1994. Ocean plankton populations as excitable media. Bull. Math. Biol. 56: 981998.). This leads us to analyse the structure of the stability or Jacobian matrix J (also known as the ‘community matrix’ (May 1974May R.M. 1974. Stability and Complexity in Model Ecosystems. Princeton University Press. 263 pp.)), composed of the partial derivatives of the original model equations (e.g. ∂(dP_{i} /dt)/∂P_{i}), ∂(dZ/dt)/∂P_{i}, etc.). In our case we have a 3×3 stability matrix given by
$$J=\left(\begin{array}{ccc}{\widehat{P}}_{11}& {\widehat{P}}_{12}& {\widehat{P}}_{1z}\\ {\widehat{P}}_{21}& {\widehat{P}}_{22}& {\widehat{P}}_{2z}\\ {\widehat{Z}}_{1}& {\widehat{Z}}_{2}& {\widehat{Z}}_{z}\end{array}\right)$$ 
(9) 
where
$${\widehat{P}}_{11}=\frac{d{B}_{1}}{d{P}_{1}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial {P}_{1}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\widehat{P}}_{22}=\frac{d{B}_{2}}{d{P}_{2}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{2}}$$ 
(10) 
${\widehat{P}}_{12}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial {P}_{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\widehat{P}}_{1z}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{1}}{\partial Z}$ 
(11) 
$${\widehat{P}}_{21}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial {P}_{1}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\widehat{P}}_{2z}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {D}_{2}}{\partial Z}$$ 
(12) 
$$\begin{array}{c}{\widehat{Z}}_{1}=\frac{\partial {D}_{1}}{\partial {P}_{1}}+\frac{\partial {D}_{2}}{\partial {P}_{1}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\widehat{Z}}_{2}=\frac{\partial {D}_{1}}{\partial {P}_{2}}+\frac{\partial {D}_{2}}{\partial {P}_{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\\ {\widehat{Z}}_{z}=\frac{\partial {D}_{1}}{\partial Z}+\frac{\partial {D}_{2}}{\partial Z}\frac{\partial E}{\partial Z}\end{array}$$ 
(13) 
Applying the relationships I we obtain that $\widehat{Z}$ _{1} = $\widehat{Z}$_{2} = 0. Then, the characteristic equation of J would be given by J − λ Id = 0, where J − λ Id means that the determinant of the matrix, Id, is the identity matrix and λ the eigenvalues of the system. Thus we have
$$\leftJ\lambda Id\right=\left(\begin{array}{ccc}{\widehat{P}}_{11}\lambda & {\widehat{P}}_{12}& {\widehat{P}}_{1z}\\ {\widehat{P}}_{21}& {\widehat{P}}_{22}\lambda & {\widehat{P}}_{2z}\\ 0& 0& {\widehat{Z}}_{z}\lambda \end{array}\right)=0$$ 
(14) 
which leads to a cubic polynomial equation in λ
λ_{3}−λ_{2}($\widehat{P}$_{11} + $\widehat{P}$_{22} + $\widehat{Z}$_{z})−λ($\widehat{P}$_{21} $\widehat{P}$_{12} − $\widehat{P}$_{11} $\widehat{Z}$_{z} − $\widehat{P}$_{22} $\widehat{Z}$_{z} −$\widehat{P}$_{11} $\widehat{P}$_{22}) > – $\widehat{P}$_{11} $\widehat{P}$_{22} $\widehat{Z}$_{z} + $\widehat{P}$_{21} $\widehat{P}$_{12} $\widehat{Z}$_{z} = 0 
(15) 
An equilibrium solution (P_{1}^{*}, P_{2}^{*}, Z ^{*}) will be stable if it satisfies the RouthHurwitz conditions (RH) (Murray 1989May R.M. 1974. Stability and Complexity in Model Ecosystems. Princeton University Press. 263 pp.). The first necessary RH condition for the stability of the system imposes that all the coefficients in the characteristic equation (15) must be positive. By contrast, if at least one of the coefficients is negative, the solution is unstable. We are interested in the case of prebloom conditions, where phytoplankton populations are typically assumed to have a low density, so it is reasonable to consider them to be positive (P_{1}, P_{2}, Z)>(0, 0, 0), but remaining near the origin. A bloom will develop when one population explodes and moves away from the origin.
Therefore, it is interesting to know which constraints determine the instability of the system in this region. According to the RH conditions, this will happen when
− $\widehat{P}$ _{11} $\widehat{P}$_{22} $\widehat{Z}$_{z} + $\widehat{P}$_{21} $\widehat{P}$_{12} $\widehat{Z}$_{z} < 0 
(16) 
or
( $\widehat{P}$ _{21} $\widehat{P}$ _{12} − $\widehat{P}$ _{11} $\widehat{Z}$ _{z} − $\widehat{P}$ _{22} $\widehat{Z}$ _{z} − $\widehat{P}$ _{11} $\widehat{P}$ _{22}) > 0 
(17) 
or
( $\widehat{P}$ _{11} + $\widehat{P}$ _{22} + $\widehat{Z}$ _{z}) > 0 
(18) 
The term involved in the last inequality (18) will always be positive because the three populations are positive and should increase ($\widehat{P}$_{11}>0, $\widehat{P}$_{22}>0, $\widehat{Z}$_{z}>0) for low values. Recall that under prebloom situations we have initially low populations levels with favourable environmental conditions (such as a sufficient nutrient supply). Therefore it is reasonable to assume that after a small perturbation both P_{1} and P_{2} will start to grow and, as a response, the grazer will tend to increase (also in agreement with the diagnostic of Gentleman et al. (2003)Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875. for functions of Class 1 and 2). As RH criteria require that all the coefficients in the polynomial must be positive, if (18) holds then at least one of the coefficients will be negative, implying that a system following the relationships IIII is unstable.
RESULTSTop
As explained above, relationships IIII are independent of the detailed functions or parameters chosen. Next, we will apply the conclusions of the previous sections to models including selected grazing functions, in order to ascertain which combination of parameters would produce blooms of a dominant species with a lower intrinsic growth rate than the nondominant one.
Parameter relationships for selected grazing functions
We selected an example of grazing function from Class 1 (MichaelisMenten or Disk) and another from Class 2 (Sigmoidal I) of the classification of Gentleman et al. (2003)Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875. (see Table 1). The first expression characterizes a constant prey preference for herbivores and the second is often used to represent passive switching (see Gentleman et al. 2003Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875., Gentleman and Neuheimer 2008Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, foodlimitation and acclimation. J. Plank. Res. 30: 12151231.; for details). In the chosen MichaelisMenten and Sigmoidal expressions, mi is the maximum ingestion rate on prey i and a_{i} = t_{i} h_{i} k, where t_{i} is the attack rate on prey i, h_{i} is the corresponding handling time and k is the halfsaturation constant for all preys (Gentleman et al. 2003Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875.).
Function: D_{i} (P_{1}, P_{2})  $$\frac{\partial {D}_{i}}{\partial {P}_{i}}$$  $$\frac{\partial {D}_{i}}{\partial {P}_{j}}$$  Name  

A  $$\frac{{m}_{i}{a}_{i}{P}_{i}}{k+{a}_{1}{P}_{1}+{a}_{2}{P}_{2}}$$  >0  <0  MichaelisMenten 
B  $$\frac{{m}_{i}{a}_{i}{P}_{i}^{2}}{{k}^{2}+{a}_{i}{P}_{1}^{2}+{a}_{i}{P}_{2}^{2}}$$  >0  <0  Sigmoidal I 
Algal growth will be represented by a logistic expression: B_{i}>(P_{i} ) = µ_{i} P_{i} − η_{i} P_{i}^{2}, where µ_{i} is the intrinsic growth rate species i and η_{i} = µ_{i} /κ_{i} is the quotient between µ_{i} and the carrying capacity (κ_{i}) for species i. A linear mortality function E(Z) = δZ will be adopted for zooplankton. Writing it as an equations system,
$$\begin{array}{c}\frac{d{P}_{1}}{dt}={\mu}_{1}{P}_{1}{\eta}_{1}{P}_{1}^{2}{D}_{1}\left({P}_{1},{P}_{2},Z\right)\\ \frac{d{P}_{2}}{dt}={\mu}_{2}{P}_{2}{\eta}_{2}{P}_{2}^{2}{D}_{2}\left({P}_{1},{P}_{2},Z\right)\\ \frac{dZ}{dt}={D}_{1}\left({P}_{1},{P}_{2},Z\right)+{D}_{2}\left({P}_{1},{P}_{2},Z\right)\delta Z\end{array}$$ 
(19) 
The initial values and parameters chosen here are derived from previous works (Solé et al. 2006aSolé J., Estrada M., GarcíaLadona E. 2006a. Biological controls of Harmful Algal Blooms: A modelling study. J. Mar. Syst. 61: 165179.,bSolé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664.). We consider (P_{1}^{0} , P_{2}^{0}, Z^{ 0} )=(0.1 mg C m^{−3}, 0.1 mg C m^{−3}, 6 mg C m^{−3}) as representative of initial low population densities near the origin. The phytoplankton saturation parameters are equal η_{1}=η_{2}=0.0001 m^{3} (mg C)^{−1} day^{−1} and the mortality parameter for zooplankton is δ=0.1 day^{−1} . The intrinsic growth rate of each species is fixed at µ_{1}=0.56 day^{−1} and µ_{2}=0.49 day^{−1}, thus giving the nondominant character to P_{1} with an intrinsic growth rate greater than that of the blooming species P_{2} . After this selection of parameters, we derive the values/ranges for m_{i}, a_{i} and k that verify, using Equations 3 and 4 and Inequalities 6, 7 and 8, the relationships I, II, III with the functions shown in Table 1. The resulting expressions are written in Table 2. The calculations were carried out by means of the Matlab© symbolic tools.
$$\frac{{{}_{}}^{}}{}$$

Table 3 shows the parameter range of both functions that satisfy the relationships IIII. These ranges have been obtained by varying one parameter at a time while keeping the rest of them fixed. Values outside the established intervals do not satisfy the relationships IIII and let P_{1} bloom or result in the same growth rate for the two prey populations.
Function: D_{i} (P_{1}, P_{2})  Range of parameter values satisfying Relationships IIII 

MichaelisMenten  m_{1}≥1, m_{2}=[0.1, 0.5] k≤25, a_{1}=[1, 10], a_{2}=[1, 15] 
Sigmoidal I  m_{1}=[0.3, 1], m_{2}=[0.1, 10] k≤35, a_{1}=[6, 20], a_{2}=[0.1, 20] 
Model simulations
In order to illustrate the population dynamics of the phytoplankton species, we carried out a numerical integration of model (19) using two sets of parameters for each grazing function (Table 4), one of them within the range given in Table 3 and triggering the dominance of P_{2} over P_{1} and the other one outside this range and leading to dominance of P_{1} over P_{2}. Figure 1 shows the time evolution of the difference P_{2} − P_{1}. If the parameters of the grazing function satisfy the relationships IIII, the net growth rate of P_{2} is higher than the net growth rate of P_{1} (Fig. 1A), while the opposite happens with the other parameter set (Fig. 1B). The two grazing functions show similar qualitative behaviour although in the Sigmoidal case the divergence between the two populations is faster. Due to its higher nonlinearity near zero, the Sigmoidal function amplifies the differences between the two populations more effectively than the MichaelisMenten. This discrepancy in the modelled phytoplankton population dynamics highlights the importance of the choice of parameters and type of predation function.
Function: D_{i} (P_{1}, P_{2})  Satisfying Rel. IIII  Not satisfying Rel. IIII 

MichaelisMenten  m_{1}=1, m_{2}=0.5 k=15, a_{1}=1, a_{2}=1  m_{1}=1, m_{2}=1 k=15, a_{1}=1, a_{2}=2 
Sigmoidal I  m_{1}=1, m_{2}=0.1 k=15, a_{1}=20, a2 =1  m_{1}=1, m_{2}=1 k=15, a_{1}=1, a_{2}=1 
In order to test the model setup, we applied the same parameters and initial conditions, setting one of the phytoplankton species to zero as an initial condition (Fig. 1C and D); as we can see, the other species grows unbounded in both cases.
To examine whether these constraints, obtained using a simple threespecies model, can be applied to a more complex ecosystem model under realistic environmental conditions, we considered simulations with the European Regional Seas Ecosystem Model (ERSEM). This model (Baretta et al. 1995BarettaBekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363379., Baretta and BarettaBekker 1997Baretta J., BarettaBekker J. 1997. Special issue: European regional seas ecosystem model II. J. Sea Res. 38: 169413., BarettaBekker et al. 1995BarettaBekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363379., Radach and Pätsch 1997Radach G., Pätsch J. 1997. Climatological annual cycles of nutrients and chlorophyll in the North Sea. J. Sea Res. 38: 231248., Pätsch and Radach 1997Pätsch J., Radach G. 1997. Longterm simulation of the eutrophication of the North Sea: temporal development of nutrients, chlorophyll and primary production in comparison to observations. J. Sea Res. 38: 275310.) is composed of interlinked modules describing biological and chemical processes within coupled water column and benthic systems. Here, however, we have used a reduced version including only the pelagic submodel. The explicit physical forcing of ERSEM includes irradiance, water temperature and wind mixing of the water column. The model characterizes planktonic functional groups in terms of physiological processes (ingestion, respiration, excretion, egestion, etc.) and uses a food web matrix to define trophic interactions. The functional groups of primary producers in the original model include picoplankton, diatoms and autotrophic flagellates (Ebenhöh et al. 1997Ebenhöh W., BarettaBekker J., Baretta J. 1997. The primary production module in the marine ecosystem model ERSEM II, with emphasis on the light forcing. J. Sea Res. 38: 173193.).
To test our conditions, we added another functional group of autotrophic flagellates (hereafter referred to as Flagellates 2 or P_{2}) formed by flagellates with a lower intrinsic growth rate but otherwise with identical characteristics to the autotrophic flagellates defined originally in the model (Flagellates 1 or P_{1}). The herbivore groups include micro and mesozooplankton and their phytoplankton grazing function is of a MichaelisMenten type. The system was forced with a time series of irradiance, water temperature, wind intensity and nutrient concentrations measured during one year in Barcelona Harbour (see details in Solé et al. 2006bSolé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664.). Irradiance and wind conditions were measured by the Spanish National Meteorological service at sampling intervals of 6 h and water temperature and nutrient concentrations were determined twice per month. All data were interpolated to obtain a sampling rate of 1 h, which is the resolution used for the time integration in the model. The parameters were set as in Solé et al. (2006b)Solé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664., with intrinsic growth rates of µ_{1} =2.02 day^{−1} for P_{1} and µ_{2} =1.77 day^{−1} for P_{2}. The initial concentrations of P_{2} (0.5 mg C m^{−3}) were assumed to be much smaller than those of P_{1}, (10 mg C m^{−3}) and the model was run over two years. The first oneyear period covered the spinup and initial numerical transients. Two oneyear runs were made, one with a set of parameters that, between July and October, did not satisfy the conditions of P_{2} dominance and one with a different set of parameters that fulfilled these conditions during the whole year (Table 4). It must be noted that in these simulations, several parameters, such as the intrinsic growth rates of the microalgae, are functions of temperature (or other variables) and may vary with time. Thus, in principle, the conditions derived in Table 4 should be recalculated for each time step of the model. However, in a first approach, a parameter can be considered as constant for a time interval if its variability during this period is slower than the characteristic time scale of the system (Vichi et al. 2007Vichi M., Pinardi N., Masina S. 2007. A generalized model of pelagic biogeochemistry for the global ocean ecosystem. Part I: Theory. J. Mar. Syst. 64: 89109.), which in this case we consider to be given by the inverse of the phytoplankton net growth rates (typically a few days).
The results of the simulations are shown in Figure 2. Diatoms (in green) have a bloom in winter with higher biomass than flagellates (note the different scales of the figure). Flagellates 1 (which have higher intrinsic growth rate than Flagellates 2) show a winter peak, which is typical of the seasonal pattern in the Mediterranean, while both Flagellates 1 and 2 exhibit several peaks along summer and fall. When relationships IIII are satisfied (Fig. 2A) during a time period (marked with continuous horizontal black line on top of the curves), Flagellates 2 produce larger blooms than Flagellates 1. When relationships IIII are not satisfied, the abundance of Flagellates 1 is always greater (Fig. 2B).
DISCUSSIONTop
We have used a simple threeequation system (two competing microalgal species and a zooplankton grazer) to deduce a set of relationships that the parameters of the feeding response functions should satisfy to allow a species with a lower intrinsic growth rate to become dominant over its competitor. This type of situation may apply to algal blooms in coastal and estuarine areas with high nutrient levels. Such events occur often. In this context, changes in selective grazing pressure can be a key element triggering the dominance of one particular species over the others. The ecological implications of differences in edibility of phytoplankton preys exposed to grazing have been explored by several models (Leibold 1989Leibold M.A. 1989. Resource edibility and the effects of predators and productivity on the outcome of trophic interactions. Am. Nat. 134: 922949., Kretzschmar et al. 1993Kretzschmar M., Nisbet R., MacCauley E. 1993. A predatorprey model for zooplankton grazing on competing algal populations. Theor. Popul. Biol. 44: 3266., Armstrong 1994Armstrong R.A. 1994. Grazing limitation and nutrient limitation in marine ecosystems: Steady state solutions of an ecosystem model with multiple food chains. Limnol. Oceanogr. 39(3): 597608., 2003Armstrong R.A. 2003. A hybrid spectral representation of phytoplankton growth and zooplancton response: The “control rod” model of plankton interaction. Deep Sea Res. Part II 50: 28952916., Holt et al. 1994Holt R., Grover J., Tilman D. 1994. Simple rules for interspecific dominance in systems with exploitative and apparent competition. Am. Nat. 144(5): 741771., Grover 1995Grover J.P. 1995. Competition, herbivory, and enrichment: nutrientbased models for edible and inedible plants. Am. Nat. 145(5): 746774.) and experiments (Bell 2002Bell T. 2002. The ecological consequences of unpalatable prey: phytoplankton response to nutrient and predator additions. Oikos 99(1): 5968.). In general, these theoretical approaches focus on the equilibrium and stability of the coexistence of predatorprey components.
Here, we centre our attention on the transient dynamics of the system. Algal blooms often appear as perturbations that move the system out of the equilibrium state in such a way that one of the species becomes dominant. Our stability analysis has shown that a system satisfying the relationships IIII must be unstable, highlighting the fact that an algal bloom may also be understood as an destabilization of the system, as assumed in models based on the theory of excitable media, which relate bloom development to the attainment of a critical value of the net growth rate of the blooming species (Truscott 1995Truscott J. 1995. Environmental forcing of simple plankton models. J. Plankton Res. 17: 22072232., Pitchford and Brindley 1999Pitchford J., Brindley J. 1999. Iron limitation, grazing pressure and oceanic high nutrientlow chlorophyll (HNLC) regions. J. Plankton Res. 21: 525547.).
Our model assumed that the initial states are characterized by low levels of phytoplankton populations, as typical of preblooms situations, and explored the interesting case in which one species with lower intrinsic growth rate than another may, however, be able to bloom due to an enhanced resistance to grazers. Given that the choice of grazing functional responses and the value of their parameters is a key point in a model construction, we applied our analysis to specific grazing functions (Teramoto et al. 1979Teramoto E., Kawasaki K., Shigesada N. 1979. Switching effect of predation on competitive prey species. J. Theor. Biol. 79: 303315., Gentleman et al. 2003Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875., Gentleman and Neuheimer 2008Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, foodlimitation and acclimation. J. Plank. Res. 30: 12151231., Anderson et al. 2010Anderson T.R., Gentleman W.C., Sinha B. 2010. Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model. Progr. Oceanogr. 87(14): 201213.). Our choice of MichaelisMenten and Sigmoidal I as examples of grazing functions is justified by the different types of behaviour (noswitching and passive switching) shown by the two functions. As shown by numerical simulations, the higher nonlinearity of the Sigmoidal with respect to the MichaelisMenten function at low phytoplankton concentrations induces a faster divergence between the two phytoplankton populations. As explained in the introduction, ecological tradeoffs between intrinsic growth rate and presence of grazing defence mechanisms have been demonstrated in association with characteristics such as nutrient composition of the prey (Yoshida et al. 2004Yoshida T., Hairston N.H., Ellner S. 2004. Evolutionary tradeoff between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris. Proc. R. Soc. London 271: 19471953., Flynn et al. 1996Flynn K.J., Davidson K., Cunningham A. 1996. Prey selection and rejection by a microflagellate; implications for the study and operation of microbial food webs. J. Exp. Mar. Biol. Ecol. 196: 357372.) and the presence of morphological defences (van Donk 1997van Donk E. 1997. Defenses in phytoplankton against grazing induced by nutrient limitation, UVB stress and infochemicals. Aq. Ecol. 31: 5358.). Toxin production by HABforming species has been shown to provide ecological advantages (Guisande et al. 2002Guisande C., Frangópulos M., Maneiro I., et al. 2002. Ecological advantages of toxin production by the dinoflagellate Alexandrium minutum under phosphorus limitation. Mar. Ecol. Prog. Ser. 225: 169176.) and it is conceivable that the production of certain chemical defences could represent a cost in terms of a lower growth rate. However, we have no data to support this possibility.
Our results are based on a simplified system, with one grazer and two preys interacting through grazing functions. A more rigorous analysis should take into account more complex models with nutrient dynamics and competition, but this would complicate the mathematical treatment. In this context, we tested the wider applicability of the relationships used in our threeequation model by employing them in ERSEM, a multispecies ecological model. The ERSEM simulations showed that when the parameters of the functional response functions fulfilled the constraints deduced for the simplified model (relationships IIII), the microalga with the lower intrinsic growth rate could grow faster than its competitor. The fact that the relationships deduced in a threedimensional system held in a higher dimensional model does not constitute a validation of them, but suggests that under the conditions applied in our ERSEM simulations, the ecological interactions between P_{1}, P_{2} and the grazers determined the behaviour of these populations (Schaffer 1981Schaffer W.M. 1981. Ecological Abstraction: The Consequences of Reduced Dimensionality in Ecological Models. Ecol. Monogr. 51: 383401.). This allows us to hypothesize that, under favourable conditions, in practice, the system dimensionality is reduced to one with three variables, as in the system we have explored here.
ACKNOWLEDGEMENTSTop
We acknowledge the support from the BIOHAB (Biological Control of Harmful Algal Blooms in European Coastal waters) project funded by the European Commission Programme on Environment and Sustainable Development of Marine Ecosystems (ESD, FPV, Key Action 3) under contract EVK3CT199900015. JS acknowledges a CSIC JAEDoc contract cofunded by the FSE. We wish to thank the two anonymous referees and the editor for their helpful comments on the manuscript.
REFERENCESTop
Anderson T.R., Gentleman W.C., Sinha B. 2010. Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model. Progr. Oceanogr. 87(14): 201213.
http://dx.doi.org/10.1016/j.pocean.2010.06.003
Armstrong R.A. 1994. Grazing limitation and nutrient limitation in marine ecosystems: Steady state solutions of an ecosystem model with multiple food chains. Limnol. Oceanogr. 39(3): 597608.
http://dx.doi.org/10.4319/lo.1994.39.3.0597
Armstrong R.A. 2003. A hybrid spectral representation of phytoplankton growth and zooplancton response: The “control rod” model of plankton interaction. Deep Sea Res. Part II 50: 28952916.
http://dx.doi.org/10.1016/j.dsr2.2003.07.003
Baretta J., BarettaBekker J. 1997. Special issue: European regional seas ecosystem model II. J. Sea Res. 38: 169413.
http://dx.doi.org/10.1016/S13851101(97)000543
Baretta J.W., Ebenhöh W., Ruardij P. 1995. The European regional seas ecosystem model: a complex marine ecosystem model. Neth. J. Sea Res. 33: 233246.
http://dx.doi.org/10.1016/00777579(95)900470
BarettaBekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363379.
http://dx.doi.org/10.1016/00777579(95)900535
Bell T. 2002. The ecological consequences of unpalatable prey: phytoplankton response to nutrient and predator additions. Oikos 99(1): 5968.
http://dx.doi.org/10.1034/j.16000706.2002.990106.x
Buskey E., Hyatt C. 1995. Effects of the Texas (USA) ‘brown tide’ alga on planktonic grazers. Mar. Ecol. Prog. Ser. 126: 285292.
http://dx.doi.org/10.3354/meps126285
Cropp R., Norbury J. 2013. Modelling plankton ecosystems and the Library of Lotka. In: Blackford J., Allen I., et al., Advances in Marine Ecosystem Modelling Research (III). J. Mar. Syst. 125: 313.
http://dx.doi.org/10.1016/j.jmarsys.2012.08.005
Ebenhöh W., BarettaBekker J., Baretta J. 1997. The primary production module in the marine ecosystem model ERSEM II, with emphasis on the light forcing. J. Sea Res. 38: 173193.
http://dx.doi.org/10.1016/S13851101(97)000439
Flynn K.J. 2002. Toxin production in migrating dinoagellates: a modelling study of PSP producing alexandrium. Harmful Algae, 1: 147155.
http://dx.doi.org/10.1016/S15689883(02)000288
Flynn K.J. 2010. Do external resource ratios matter? Implications for modelling eutrophication events and controlling harmful algal blooms. J. Mar. Syst. 83: 170180.
http://dx.doi.org/10.1016/j.jmarsys.2010.04.007
Flynn K.J., Iringoien X. 2009. Why aldehydeinduced insidious effects cannot be considered as a diatom defence mechanism against copepods. Mar. Ecol. Prog. Ser. 377: 7989.
http://dx.doi.org/10.3354/meps07865
Flynn K.J., Davidson K., Cunningham A. 1996. Prey selection and rejection by a microflagellate; implications for the study and operation of microbial food webs. J. Exp. Mar. Biol. Ecol. 196: 357372.
http://dx.doi.org/10.1016/00220981(95)001409
Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, foodlimitation and acclimation. J. Plank. Res. 30: 12151231.
http://dx.doi.org/10.1093/plankt/fbn078
Gentleman W., Leising A., Frost B., et al. 2003. Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics. Deep Sea Res. Part II 50: 28472875.
http://dx.doi.org/10.1016/j.dsr2.2003.07.001
Grover J.P. 1995. Competition, herbivory, and enrichment: nutrientbased models for edible and inedible plants. Am. Nat. 145(5): 746774.
http://dx.doi.org/10.1086/285766
Guisande C., Frangópulos M., Maneiro I., et al. 2002. Ecological advantages of toxin production by the dinoflagellate Alexandrium minutum under phosphorus limitation. Mar. Ecol. Prog. Ser. 225: 169176.
http://dx.doi.org/10.3354/meps225169
Holling C. 1959. Some characteristics of simple types of predation and parasitism. Can. Entom. 91: 385398.
http://dx.doi.org/10.4039/Ent913857
Holt R., Grover J., Tilman D. 1994. Simple rules for interspecific dominance in systems with exploitative and apparent competition. Am. Nat. 144(5): 741771.
http://dx.doi.org/10.1086/285705
Irigoien X., Flynn K.J., Harris R.P. 2005. Phytoplankton blooms: a loophole in microzooplankton grazing impact? J. Plankton Res. 27: 313321.
http://dx.doi.org/10.1093/plankt/fbi011
Jessup C.M., Bohannan B.J.M. 2008. The shape of an ecological tradeoff varies with environment. Ecol. Let. 11: 947959.
http://dx.doi.org/10.1111/j.14610248.2008.01205.x
Kierstead H., Slobodkin L. 1953. The size of water masses containing plankton blooms. J. Mar. Res. 12: 141147.
Kretzschmar M., Nisbet R., MacCauley E. 1993. A predatorprey model for zooplankton grazing on competing algal populations. Theor. Popul. Biol. 44: 3266.
http://dx.doi.org/10.1006/tpbi.1993.1017
Leibold M.A. 1989. Resource edibility and the effects of predators and productivity on the outcome of trophic interactions. Am. Nat. 134: 922949.
http://dx.doi.org/10.1086/285022
Leibold M.A. 1996. A graphical model of keystone predators in food webs: trophic regulation of abundance, incidence, and diversity patterns in communities. Am. Nat. 147: 784812.
http://dx.doi.org/10.1086/285879
Leising A., Horner R., Pierson J., et al. 2005. The balance between microzooplankton grazing and phytoplankton growth in a highly productive estuarine fjord. Progr. Oceanogr. 67: 366383.
http://dx.doi.org/10.1016/j.pocean.2005.09.007
Litchman E., Klausmeier C. 2008. Traitbased community ecology of phytoplankton. Annu. Rev. Ecol. Evol. Syst. 39: 615639.
http://dx.doi.org/10.1146/annurev.ecolsys.39.110707.173549
Maestrini S., Granéli E. 1991. Environmental conditions and ecophysiological mechanisms which led to the 1988 Chrysochromulina polylepis bloom: an hypothesis. Ocean. Acta 14: 397413.
Margalef R., Estrada M., Blasco D. 1979. Functional morphology of organisms involved in red tides, as adapted do cecaying turbulence. In: Taylor D.L. and Sliger H.H. (eds), Toxic dinoagellate blooms. Elsevier, North Holland, pp. 315320.
May R.M. 1974. Stability and Complexity in Model Ecosystems. Princeton University Press. 263 pp.
Mitra A., Flynn K. 2005. Predatorprey interactions: is ‘ecological stoichiometry’ sufficient when good food goes bad? J. Plankton Res. 27: 393399.
http://dx.doi.org/10.1093/plankt/fbi022
Mitra A., Flynn K. 2006. Promotion of harmful algal blooms by zooplankton predatory activity. Biol. Lett. 2: 194197.
http://dx.doi.org/10.1098/rsbl.2006.0447
Murray J. 1989. Mathematical Biology. Springer, Berlin.
http://dx.doi.org/10.1007/9783662085394
Pätsch J., Radach G. 1997. Longterm simulation of the eutrophication of the North Sea: temporal development of nutrients, chlorophyll and primary production in comparison to observations. J. Sea Res. 38: 275310.
http://dx.doi.org/10.1016/S13851101(97)000518
Pitchford J., Brindley J. 1999. Iron limitation, grazing pressure and oceanic high nutrientlow chlorophyll (HNLC) regions. J. Plankton Res. 21: 525547.
http://dx.doi.org/10.1093/plankt/21.3.525
Radach G., Pätsch J. 1997. Climatological annual cycles of nutrients and chlorophyll in the North Sea. J. Sea Res. 38: 231248.
http://dx.doi.org/10.1016/S13851101(97)000488
Sailley S., Vogt M., Doney S., et al. 2013. Comparing food web structures and dynamics across a suite of global marine ecosystem models. Ecol. Modell. 261262: 4357.
http://dx.doi.org/10.1016/j.ecolmodel.2013.04.006
Schaffer W.M. 1981. Ecological Abstraction: The Consequences of Reduced Dimensionality in Ecological Models. Ecol. Monogr. 51: 383401.
Smayda T.J., Reynolds C.S. 2001. Community assembly in marine phytoplankton: application of recent models to harmful dinoagellate blooms. J. Plankton Res. 23(5): 447461.
http://dx.doi.org/10.1093/plankt/23.5.447
Solé J., Estrada M., GarcíaLadona E. 2006a. Biological controls of Harmful Algal Blooms: A modelling study. J. Mar. Syst. 61: 165179.
http://dx.doi.org/10.1016/j.jmarsys.2005.06.004
Solé J., GarcíaLadona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 4664.
http://dx.doi.org/10.1016/j.jmarsys.2006.04.002
Steele J. 1974. The structure of marine ecosystems. Harvard University Press.
http://dx.doi.org/10.4159/harvard.9780674592513
Teramoto E., Kawasaki K., Shigesada N. 1979. Switching effect of predation on competitive prey species. J. Theor. Biol. 79: 303315.
http://dx.doi.org/10.1016/00225193(79)903485
Tillmann U., Hesse K., Colijn F. 2000. Planktonic primary production in the German Wadden Sea. J. Plankton Res. 22: 12531276.
http://dx.doi.org/10.1093/plankt/22.7.1253
Touzet N., Franco J., Raine R. 2007. Influence of inorganic nutrition on growth and PSP toxin production of Alexandrium minutum (Dinophyceae) from Cork Harbour, Ireland. Toxicon 50: 106119.
http://dx.doi.org/10.1016/j.toxicon.2007.03.001
Truscott J. 1995. Environmental forcing of simple plankton models. J. Plankton Res. 17: 22072232.
http://dx.doi.org/10.1093/plankt/17.12.2207
Truscott J., Brindley J. 1994. Ocean plankton populations as excitable media. Bull. Math. Biol. 56: 981998.
http://dx.doi.org/10.1007/BF02458277
van Donk E. 1997. Defenses in phytoplankton against grazing induced by nutrient limitation, UVB stress and infochemicals. Aq. Ecol. 31: 5358.
http://dx.doi.org/10.1023/A:1009951622185
Vichi M., Pinardi N., Masina S. 2007. A generalized model of pelagic biogeochemistry for the global ocean ecosystem. Part I: Theory. J. Mar. Syst. 64: 89109.
http://dx.doi.org/10.1016/j.jmarsys.2006.03.006
Wyatt T., Horwood J. 1973. Model which generates red tides. Nature 244: 238240.
http://dx.doi.org/10.1038/244238a0
Yoshida T., Hairston N.H., Ellner S. 2004. Evolutionary tradeoff between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris. Proc. R. Soc. London 271: 19471953.
http://dx.doi.org/10.1098/rspb.2004.2818