SUMMARY

SCIENTIA MARINA

Scientia Marina

Sci Mar

0214-8358

Consejo Superior de Investigaciones Científicas

sm4271

10.3989/scimar.04271.18A

Parameter constraints of grazing response functions. Implications for phytoplankton bloom initiation

Limitaciones de los parámetros de las funciones de predación. Implicaciones para la iniciación de las proliferaciones de fitoplancton

Parameter constraints for bloom initiation

Pelegrí

Josep L.

Special Issue Editor Vaqué

Dolors

Special Issue Editor

Solé

Jordi

García-Ladona

Emilio

Piera

Jaume

Estrada

Marta

Institut de Ciències del Mar, CSIC, 08003 Barcelona, Catalunya, Spain.

e-mail: jsole@icm.cat

30 09 2016

2016

80S1 Suppl. 1

Planet Ocean

129 137

10.3989/scimar.04271.18A

14 5 2015 15 9 2015 30 9 2016

2016

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-by) Spain 3.0 License.

SUMMARY

Phytoplankton blooms are events of production and accumulation of phytoplankton biomass that influence ecosystem dynamics and may also have effects on socio-economic activities. Among the biological factors that affect bloom dynamics, prey selection by zooplankton may play an important role. Here we consider the initial state of development of an algal bloom and analyse how a reduced grazing pressure can allow an algal species with a lower intrinsic growth rate than a competitor to become dominant. We use a simple model with two microalgal species and one zooplankton grazer to derive general relationships between phytoplankton growth and zooplankton grazing. These relationships are applied to two common grazing response functions in order to deduce the mathematical constraints that the parameters of these functions must obey to allow the dominance of the lower growth rate competitor. To assess the usefulness of the deduced relationships in a more general framework, the results are applied in the context of a multispecies ecosystem model (ERSEM).

RESUMEN

Las proliferaciones de fitoplancton son eventos de producción y acumulación de biomasa de fitoplancton que tienen una fuerte influencia tanto en la dinámica del ecosistema como en actividades socioeconómicas. Entre los factores biológicos que afectan la dinámica de las proliferaciones fitoplanctónicas, la selección de presas por el zooplancton puede jugar un papel importante. En este trabajo consideramos el estado inicial de desarrollo de una proliferación algal considerando que la presión de la predación puede permitir a una especie de algas con una menor tasa de crecimiento intrínseco que un competidor, ser dominante. Utilizamos un modelo sencillo con dos especies de microalgas y un herbívoro (zooplancton) para deducir relaciones generales entre el crecimiento del fitoplancton y la predación. Estas relaciones se aplican a dos funciones de respuesta de predación para deducir las limitaciones matemáticas que los parámetros de estas funciones deben obedecer. Para evaluar la utilidad de las relaciones deducidas en un marco más general, los resultados se aplican a un modelo de múltiples especies del ecosistema marino (ERSEM).

KEYWORDS algal blooms plankton prey selection grazing functions multispecies model mathematical constraints

PALABRAS CLAVE proliferaciones algales plancton selección de la presa funciones de predación modelo multiespecie limitaciones matemáticas

INTRODUCTION

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 socio-economic wellbeing of the area in which the bloom occurs. Basic physico-chemical conditions for bloom development include nutrient availability and low physical dispersion losses (Kierstead and Slobodkin 1953, Wyatt and Horwood 1973, Margalef et al. 1979, Smayda and Reynolds 2001). 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 1974, Holt et al. 1994, Leibold 1996, Armstrong 2003). Prey selection may allow or inhibit the proliferation of certain phytoplankton species, as has been shown for high nutrient-low chlorophyll (HNLC) regions (Leising et al. 2005). Some harmful algal blooms (HABs) have been attributed to grazing impairment due to toxin production. Well-known examples include the Chrysochromulina polylepis red tide studied by Maestrini and Granéli (1991) 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 1995). 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. 2005).

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 trade-offs can have important impacts on food web dynamics and have been the subject of numerous studies (Leibold 1989, Tillmann et al. 2000, Jessup and Bohannan 2008). Trade-offs 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. 2004, Litchman and Klausmeier 2008). Another example would be the metabolic costs associated with the production of morphological defences (van Donk 1997). The importance of ’prey stoichiometry’ in predator-prey interactions has been well documented (e.g. Flynn et al. 1996, Mitra and Flynn 2005, 2006, Flynn 2010). Nutrient stress may also interact with growth rate and toxin production (Flynn 2002, Flynn and Iringoien 2009, Touzet et al. 2007). The role of grazing-deterrent toxins in allowing a toxic algal species to develop a bloom in competition with a non-defended competitor with higher growth rate was explored in a general framework by Solé et al. (2006b), by means of a simple predator-prey model. These authors used a logistic function for algal growth and a Holling II–type grazing response (Holling 1959) 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 non-defended 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), Sailley et al. (2013) and Cropp and Norbury (2013), 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) 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 pre-bloom 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. 2006a,b). 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 three-species 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 METHODS

Because our main interest is pre-bloom conditions, we will focus on the short time response starting from a near-equilibrium 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 ‘non-dominant’ species) with a higher intrinsic growth rate (of course, the case is trivial if the ‘dominant’ species has a higher growth rate than the ‘non-dominant’ one).

General conditions on feeding response functions

We propose a simplified ecosystem model that includes two microalgal populations, a non-dominant (P₁) and a dominant one (P₂), 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):

\begin{matrix} \frac{d P_{1}}{d t} = B_{1} (P_{1}) - D_{1} (P_{1}, P_{2}, Z) \\ \frac{d P_{2}}{d t} = B_{2} (P_{2}) - D_{2} (P_{1}, P_{2}, Z) \\ \frac{d Z}{d t} = D_{1} (P_{1}, P_{2}, Z) + D_{2} (P_{1}, P_{2}, Z) - E (Z) \end{matrix}

(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 pre-bloom stage the nutrients do not have a limiting effect on algal growth.

In principle, all functions and parameters should be time-dependent following physiological or environmental changes. However, as our study will address the short-term response of the system, starting from pre-bloom conditions with low biomass concentration, we will assume that the system parameters remain constant at these pre-bloom conditions, before interspecific competition or environmental limitations become significant.

We consider that the pre-bloom system is near an equilibrium point (P₁^*, P₂^*, Z^*), in which source and loss rates are balanced so populations do not vary over time:

(\frac{d P_{1}}{d t} = \frac{d P_{2}}{d t} = \frac{d Z_{1}}{d t} = 0)

(2)

This state satisfies

B₁ (P₁^*) = D₁ (P₁^*, P₂^*, Z ^*)

(2a)

B₂ (P₂^*) = D2 (P₁^*, P₂^*, Z ^*)

(2b)

E(Z^*) = D₁ (P₁^*, P₂^*, Z ^*) + D₂(P₁^*, P₂^*, 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₁ 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 or \frac{\partial D_{2}}{\partial P_{1}} and \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₁ and D₂ do not depend on P₁, implying that grazers control only the growth of P₂, leaving P₁ unbounded by grazing pressure. However we have assumed that both grazing functions depend on P₁ and P₁; therefore,

\frac{\partial D_{1}}{\partial P_{1}} and \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₂, 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}} and \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).

– 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}}{d t} ≃ 0, \frac{P_{2}}{d t} ≃ 0

A necessary condition for the dominant species, in our case P₂, to bloom is to have a greater net growth rate than P₁

\frac{P_{1}}{d t} < \frac{P_{2}}{d t} ≃ 0

According to the system (1), we have

B₁ − D₁ < B₂ − D₂,

(5)

Taking the derivatives (on P₁ and P₂) of this expression and using I we can obtain the relationships:

\frac{d B_{2}}{d P_{2}} >2 \frac{\partial D_{2}}{\partial P_{2}}

(6)

\frac{d B_{1}}{d P_{1}} < 2 \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}} > \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₁ must grow faster in the P₁ direction (the P₁ axis in the phase space) than the grazing function D₂ in the P₂ 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 1989, Truscott and Brindley 1994). This leads us to analyse the structure of the stability or Jacobian matrix J (also known as the ‘community matrix’ (May 1974)), 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 = (\begin{matrix} {\hat{P}}_{11} & {\hat{P}}_{12} & {\hat{P}}_{1 z} \\ {\hat{P}}_{21} & {\hat{P}}_{22} & {\hat{P}}_{2 z} \\ {\hat{Z}}_{1} & {\hat{Z}}_{2} & {\hat{Z}}_{z} \end{matrix})

(9)

where

{\hat{P}}_{11} = \frac{d B_{1}}{d P_{1}} - \frac{\partial D_{1}}{\partial P_{1}}, {\hat{P}}_{22} = \frac{d B_{2}}{d P_{2}} - \frac{\partial D_{2}}{\partial P_{2}}

(10)

{\hat{P}}_{12} = - \frac{\partial D_{1}}{\partial P_{2}}, {\hat{P}}_{1 z} = - \frac{\partial D_{1}}{\partial Z}

(11)

{\hat{P}}_{21} = - \frac{\partial D_{2}}{\partial P_{1}}, {\hat{P}}_{2 z} = - \frac{\partial D_{2}}{\partial Z}

(12)

\begin{matrix} {\hat{Z}}_{1} = \frac{\partial D_{1}}{\partial P_{1}} + \frac{\partial D_{2}}{\partial P_{1}}, {\hat{Z}}_{2} = \frac{\partial D_{1}}{\partial P_{2}} + \frac{\partial D_{2}}{\partial P_{2}}, \\ {\hat{Z}}_{z} = \frac{\partial D_{1}}{\partial Z} + \frac{\partial D_{2}}{\partial Z} - \frac{\partial E}{\partial Z} \end{matrix}

(13)

Applying the relationships I we obtain that $\hat{Z}$ ₁ = $\hat{Z}$ ₂ = 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

| J - λ I d | = (\begin{matrix} {\hat{P}}_{11} - λ & {\hat{P}}_{12} & {\hat{P}}_{1 z} \\ {\hat{P}}_{21} & {\hat{P}}_{22} - λ & {\hat{P}}_{2 z} \\ 0 & 0 & {\hat{Z}}_{z} - λ \end{matrix}) = 0

(14)

which leads to a cubic polynomial equation in λ

λ₃−λ₂(

\hat{P}

₁₁ +

\hat{P}

₂₂ +

\hat{Z}

_z)−λ(

\hat{P}

₂₁

\hat{P}

₁₂ −

\hat{P}

₁₁

\hat{Z}

_z −

\hat{P}

₂₂

\hat{Z}

_z −

\hat{P}

₁₁

\hat{P}

₂₂) > –

\hat{P}

₁₁

\hat{P}

₂₂

\hat{Z}

_z +

\hat{P}

₂₁

\hat{P}

₁₂

\hat{Z}

_z = 0

(15)

An equilibrium solution (P₁^*, P₂^*, Z ^*) will be stable if it satisfies the Routh-Hurwitz conditions (RH) (Murray 1989). 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 pre-bloom conditions, where phytoplankton populations are typically assumed to have a low density, so it is reasonable to consider them to be positive (P₁, P₂, 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

−

\hat{P}

₁₁

\hat{P}

₂₂

\hat{Z}

_z +

\hat{P}

₂₁

\hat{P}

₁₂

\hat{Z}

_z < 0

(16)

(

\hat{P}

₂₁

\hat{P}

₁₂ −

\hat{P}

₁₁

\hat{Z}

_z −

\hat{P}

₂₂

\hat{Z}

_z −

\hat{P}

₁₁

\hat{P}

₂₂) > 0

(17)

(

\hat{P}

₁₁ +

\hat{P}

₂₂ +

\hat{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 ( $\hat{P}$ ₁₁>0, $\hat{P}$ ₂₂>0, $\hat{Z}$ _z>0) for low values. Recall that under pre-bloom 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₁ and P₂ 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) 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 I-III is unstable.

RESULTS

As explained above, relationships I-III 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 non-dominant one.

Parameter relationships for selected grazing functions

We selected an example of grazing function from Class 1 (Michaelis-Menten or Disk) and another from Class 2 (Sigmoidal I) of the classification of Gentleman et al. (2003) (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. 2003, Gentleman and Neuheimer 2008; for details). In the chosen Michaelis-Menten 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 half-saturation constant for all preys (Gentleman et al. 2003).

Table 1 Selected functional responses for zooplankton feeding, based on <xref ref-type="bibr" rid="CIT16">Gentleman et al. (2003)</xref> and references therein. The parameters <italic>mi</italic> and <italic>ai</italic> are characteristic for each resource and <italic>k</italic> is the half-saturation constant of the grazing functions (see text). The columns of derivatives indicate the sign of the derivative of each grazing function with respect to <italic>P</italic>1 and <italic>P</italic>2. The columns of derivatives are as those reported in <xref ref-type="bibr" rid="CIT16">Gentleman et al. (2003)</xref>.

	Function: D_i (P₁, P₂)	$\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	Michaelis-Menten
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², 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{matrix} \frac{d P_{1}}{d t} = µ_{1} P_{1} - η_{1} P_{1}^{2} - D_{1} (P_{1}, P_{2}, Z) \\ \frac{d P_{2}}{d t} = µ_{2} P_{2} - η_{2} P_{2}^{2} - D_{2} (P_{1}, P_{2}, Z) \\ \frac{d Z}{d t} = D_{1} (P_{1}, P_{2}, Z) + D_{2} (P_{1}, P_{2}, Z) - δ Z \end{matrix}

(19)

The initial values and parameters chosen here are derived from previous works (Solé et al. 2006a,b). We consider (P₁⁰ , P₂⁰, Z⁰ )=(0.1 mg C m⁻³, 0.1 mg C m⁻³, 6 mg C m⁻³) as representative of initial low population densities near the origin. The phytoplankton saturation parameters are equal η₁=η₂=0.0001 m³ (mg C)⁻¹ day⁻¹ and the mortality parameter for zooplankton is δ=0.1 day⁻¹ . The intrinsic growth rate of each species is fixed at µ₁=0.56 day⁻¹ and µ₂=0.49 day⁻¹, thus giving the non-dominant character to P₁ with an intrinsic growth rate greater than that of the blooming species P₂ . 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.

Table 2 Conditions that parameters of the different grazing functions must satisfy to allow <italic>P</italic>2 to initiate a bloom. <italic>P</italic>10, <italic>P</italic>20 are the initial concentrations of <italic>P</italic>1 and <italic>P</italic>2. <italic>ηi</italic> = <italic>µi</italic> /<italic>κi</italic> is the quotient between the intrinsic growth rate (<italic>µi</italic>) and the carrying capacity (<italic>κi</italic>) corresponding to species <italic>i</italic>. The parameters <italic>m</italic>i and <italic>a</italic>i are characteristic of each prey and <italic>k</italic> is the half-saturation constant (see text) of the grazing functions.

Function Name	Relationship	Conditions on parameters
Michaelis-Menten	I	m₁(k + a₂P₂⁰) = m₂a₂P₂⁰
	I	m₂(k + a₁P₁⁰) = m₁a₁P₁⁰
	II	$μ_{2} + 2 η_{2} P_{2}^{0} > 2 m_{2} a_{2} \frac{k + a_{1} P_{1}^{0}}{{(k + a_{1} P_{1}^{0} + a_{2} P_{2}^{0})}^{2}}$
	II	$μ_{1} + 2 η_{1} P_{1}^{0} > 2 m_{1} a_{1} \frac{k + a_{2} P_{2}^{0}}{{(k + a_{1} P_{1}^{0} + a_{2} P_{2}^{0})}^{2}}$
	III	m₁ a₁(k + a₂P₂⁰) > m₂a₂(k + a₁P₁⁰)
Sigmoidal I	I	m₁(k₂ + a₂ (P₂⁰)²) = m₂a₂ P₂⁰
	I	m₂(k₂ + a₁ (P₁⁰)²) = m₁a₁P₁⁰
	II	$μ_{2} - η_{2} P_{2}^{0} > 4 m_{2} P_{2}^{0} a_{2} \frac{k^{2} + a_{1} {(P_{1}^{0})}^{2}}{{(k^{2} + a_{1} {(P_{1}^{0})}^{2} + a_{2} {(P_{2}^{0})}^{2})}^{2}}$
	II	$μ_{1} - η_{1} P_{1}^{0} > 4 m_{1} P_{1}^{0} a_{1} \frac{k^{2} + a_{2} {(P_{2}^{0})}^{2}}{{(k^{2} + a_{1} {(P_{1}^{0})}^{2} + a_{2} {(P_{2}^{0})}^{2})}^{2}}$
	III	m₁a₁P₁⁰ (k₂ + a₂(P₂⁰)²) > m₂ P₂⁰a₂(k₂ + a₁(P₁⁰)²)

Table 3 shows the parameter range of both functions that satisfy the relationships I-III. 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 I-III and let P₁ bloom or result in the same growth rate for the two prey populations.

Table 3 Functions and parameter value ranges of the grazing functions satisfying Relationships I-III. The ranges are established by varying one parameter at a time while keeping the others fixed. The other parameter values are <italic>µ</italic>1 = 0.56 day−1, <italic>µ</italic>2 = 0.49 day−1 <italic>η</italic>1 = <italic>η</italic>2 = 0.0001 m3 (mg C)−1 day−1 and <italic>δ</italic> = 0.1 day−1

Function: D_i (P₁, P₂)	Range of parameter values satisfying Relationships I-III
Michaelis-Menten	m₁≥1, m₂=[0.1, 0.5] k≤25, a₁=[1, 10], a₂=[1, 15]
Sigmoidal I	m₁=[0.3, 1], m₂=[0.1, 10] k≤35, a₁=[6, 20], a₂=[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₂ over P₁ and the other one outside this range and leading to dominance of P₁ over P₂. Figure 1 shows the time evolution of the difference P₂ − P₁. If the parameters of the grazing function satisfy the relationships I-III, the net growth rate of P₂ is higher than the net growth rate of P₁ (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 non-linearity near zero, the Sigmoidal function amplifies the differences between the two populations more effectively than the Michaelis-Menten. This discrepancy in the modelled phytoplankton population dynamics highlights the importance of the choice of parameters and type of predation function.

Table 4 Functions and parameter values of the grazing functions (see <xref ref-type="table" rid="T1">Table 1</xref>) used in the numerical simulations of the three-equation model (19). <italic>µ</italic>1=0.56 day−1, <italic>µ</italic>2=0.49 day−1; <italic>η</italic>1=<italic>η</italic>2=0.0001 m3 (mg C)−1 day−1 and <italic>δ</italic>=0.1 day−1.

Function: D_i (P₁, P₂)	Satisfying Rel. I-III	Not satisfying Rel. I-III
Michaelis-Menten	m₁=1, m₂=0.5 k=15, a₁=1, a₂=1	m₁=1, m₂=1 k=15, a₁=1, a₂=2
Sigmoidal I	m₁=1, m₂=0.1 k=15, a₁=20, a2 =1	m₁=1, m₂=1 k=15, a₁=1, a₂=1

Fig. 1 Difference <italic>P</italic>2−<italic>P</italic>1 as simulated with the three-equation model using the grazing functions of <xref ref-type="table" rid="T1">Table 1</xref> and the parameter values of <xref ref-type="table" rid="T4">Table 4</xref>. Initial concentrations are <italic>P</italic>1=<italic>P</italic>2=0.1 mg C m−3, Z=6 mg C m−3. A, satisfying Relationships I-III: <italic>P</italic>2 exceeds <italic>P</italic>1 (initiation of bloom by <italic>P</italic>2). B, not satisfying Relationships I-III: <italic>P</italic>1 exceeds <italic>P</italic>2 (initiation of bloom by <italic>P</italic>1). The same as the previous cases but for initial concentrations: C, <italic>P</italic>1=0 mg C m−3, <italic>P</italic>2=0.1 mg C m−3, <italic>Z</italic>=6 mg C m−3 and D, <italic>P</italic>1=0.1 mg C m−3, <italic>P</italic>2=0 mg C m−3, Z=6 mg C m−3.

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 three-species 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. 1995, Baretta and Baretta-Bekker 1997, Baretta-Bekker et al. 1995, Radach and Pätsch 1997, Pätsch and Radach 1997) 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. 1997).

To test our conditions, we added another functional group of autotrophic flagellates (hereafter referred to as Flagellates 2 or P₂) 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₁). The herbivore groups include micro- and mesozooplankton and their phytoplankton grazing function is of a Michaelis-Menten 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. 2006b). 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), with intrinsic growth rates of µ₁ =2.02 day⁻¹ for P₁ and µ₂ =1.77 day⁻¹ for P₂. The initial concentrations of P₂ (0.5 mg C m⁻³) were assumed to be much smaller than those of P₁, (10 mg C m⁻³) and the model was run over two years. The first one-year period covered the spin-up and initial numerical transients. Two one-year runs were made, one with a set of parameters that, between July and October, did not satisfy the conditions of P₂ 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. 2007), 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 I-III 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 I-III are not satisfied, the abundance of Flagellates 1 is always greater (Fig. 2B).

Fig. 2 One-year numerical simulations with the ERSEM model using the Michaelis-Menten function. Diatoms are shown in green, solid line; Flagellates 1 in black, pointed line and Flagellates 2 in red, dashed line. A, Michaelis-Menten grazing function with parameters satisfying Relationships I-III only during the July-October period (marked as a continuous black line on top of the curves); B, Michaelis-Menten grazing function with parameters not satisfying Relationships I-III. DISCUSSION

We have used a simple three-equation 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 1989, Kretzschmar et al. 1993, Armstrong 1994, 2003, Holt et al. 1994, Grover 1995) and experiments (Bell 2002). In general, these theoretical approaches focus on the equilibrium and stability of the coexistence of predator-prey 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 I-III 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 1995, Pitchford and Brindley 1999).

Our model assumed that the initial states are characterized by low levels of phytoplankton populations, as typical of pre-blooms 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. 1979, Gentleman et al. 2003, Gentleman and Neuheimer 2008, Anderson et al. 2010). Our choice of Michaelis-Menten and Sigmoidal I as examples of grazing functions is justified by the different types of behaviour (no-switching and passive switching) shown by the two functions. As shown by numerical simulations, the higher non-linearity of the Sigmoidal with respect to the Michaelis-Menten function at low phytoplankton concentrations induces a faster divergence between the two phytoplankton populations. As explained in the introduction, ecological trade-offs 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. 2004, Flynn et al. 1996) and the presence of morphological defences (van Donk 1997). Toxin production by HAB-forming species has been shown to provide ecological advantages (Guisande et al. 2002) 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 three-equation 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 I-III), the microalga with the lower intrinsic growth rate could grow faster than its competitor. The fact that the relationships deduced in a three-dimensional 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₁, P₂ and the grazers determined the behaviour of these populations (Schaffer 1981). 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.

ACKNOWLEDGEMENTS

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, FP-V, Key Action 3) under contract EVK3-CT-1999-00015. JS acknowledges a CSIC JAE-Doc contract co-funded by the FSE. We wish to thank the two anonymous referees and the editor for their helpful comments on the manuscript.

REFERENCES

Anderson

T.R.

Gentleman

W.C.

Sinha

Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model

Progr. Oceanogr. 2010 87 201 213

http://dx.doi.org/10.1016/j.pocean.2010.06.003

Armstrong

R.A.

Grazing limitation and nutrient limitation in marine ecosystems: Steady state solutions of an ecosystem model with multiple food chains

Limnol. Oceanogr. 1994 39 597 608

http://dx.doi.org/10.4319/lo.1994.39.3.0597

Armstrong

R.A.

A hybrid spectral representation of phytoplankton growth and zooplancton response: The "control rod" model of plankton interaction

Deep Sea Res. Part II 2003 50 2895 2916

http://dx.doi.org/10.1016/j.dsr2.2003.07.003

Baretta

Baretta-Bekker

Special issue: European regional seas ecosystem model II.

J. Sea Res. 1997 38 169 413

http://dx.doi.org/10.1016/S1385-1101(97)00054-3

Baretta

J.W.

Ebenhöh

Ruardij

The european regional seas ecosystem model: a complex marine ecosystem model.

Neth. J. Sea Res 1995 33 233 246

http://dx.doi.org/10.1016/0077-7579(95)90047-0

Baretta-Bekker

J.G.

Baretta

Rasmussen

The microbial food web in the european regional seas ecosystem model

Neth. J. Sea Res., 1995 33 363 379

http://dx.doi.org/10.1016/0077-7579(95)90053-5

Bell

The ecological consequences of unpalatable prey: phytoplankton response to nutrient and predator additions

Oikos 2002 99 59 68

http://dx.doi.org/10.1034/j.1600-0706.2002.990106.x

Buskey

Hyatt

Effects of the Texas(USA) 'brown tide' alga on planktonic grazers

Mar. Ecol. Prog. Ser., 1995 126 285 292

http://dx.doi.org/10.3354/meps126285

Cropp

Norbury

Advances in Marine Ecosystem Modelling Research (III).

Blackford

Allen

Modelling plankton ecosystems and the library of lotka

J. Mar. Syst. 2013 125 3 13

http://dx.doi.org/10.1016/j.jmarsys.2012.08.005

Ebenhöh

Baretta-Bekker

Baretta

The primary production module in the marine ecosystem model ERSEM II, with emphasis on the light forcing.

J. Sea Res. 1997 38 173 193

193. http://dx.doi.org/10.1016/S1385-1101(97)00043-9

Flynn

K.J.

Toxin production in migrating dinoagellates: a modelling study of psp producing alexandrium

Harmful Algae 2002 1 147 155

http://dx.doi.org/10.1016/S1568-9883(02)00028-8

Flynn

K.J.

Do external resource ratios matter? implications for modelling eutrophication events and controlling harmful algal blooms

J. Mar. Syst. 2010 83 170 180

http://dx.doi.org/10.1016/j.jmarsys.2010.04.007

Flynn

K.J.

Iringoien

Why aldehyde-induced insidious effects cannot be considered as a diatom defence mechanism against copepods.

Mar. Ecol. Prog. Ser., 2009 377 79 89

http://dx.doi.org/10.3354/meps07865

Flynn

K.J.

Davidson

Cunningham

Prey selection and rejection by a microagellate; implications for the study and operation of microbial food webs

J. Exp. Mar. Biol. Ecol. 1996 196 357 372

http://dx.doi.org/10.1016/0022-0981(95)00140-9

Gentleman

W.C.

Neuheimer

A.B.

Functional responses and ecosystem dynamics: how clearance rates explain the inuence of satiation, food-limitation and acclimation

J. Plank. Res. 2008 30 1215 1231

http://dx.doi.org/10.1093/plankt/fbn078

Gentleman

Leising

Frost

Functional responses for zooplankton feeding on multiple resources: a critical review of assumed biological dynamics

Deep Sea Res. Part II 2003 50 2847 2875

http://dx.doi.org/10.1016/j.dsr2.2003.07.001

Grover

J.P.

Competition, herbivory, and enrichment: nutrient-based models for edible and inedible plants

Am. Nat. 1995 145 746 774

http://dx.doi.org/10.1086/285766

Guisande

Frangópulos

Maneiro

Ecological advantages of toxin production by the dinoagellate Alexandrium minutum under phosphorus limitation.

Mar. Ecol. Prog. Ser. 2002 225 169 176

http://dx.doi.org/10.3354/meps225169

Holling

Some characteristics of simple types of predation and parasitism

Can. Entom. 1959 91 385 398

http://dx.doi.org/10.4039/Ent91385-7

Holt

Grover

Tilman

Simple rules for interspecific dominance in systems with explotative and apparent competition

Am. Nat. 1994 144 741 771

http://dx.doi.org/10.1086/285705

Irigoien

Flynn

K.J.

Harris

R.P.

Phytoplankton blooms: a loophole in microzooplankton grazing impact?

J. Plankton Res., 2005 27 313 321

http://dx.doi.org/10.1093/plankt/fbi011

Jessup

C.M.

Bohannan

B.J.M.

The shape of an ecological trade-off varies with environment

Ecol. Let., 2008 11 947 959

http://dx.doi.org/10.1111/j.1461-0248.2008.01205.x

Kierstead

Slobodkin

The size of water masses containing plankton blooms

J. Mar. Res., 1953 12 141 147

Kretzschmar

Nisbet

MacCauley

A predator-prey model for zooplankton grazing on competing algal populations

Theor. Popul. Biol. 1993 44 32 66

http://dx.doi.org/10.1006/tpbi.1993.1017

Leibold

M.A.

Resource edibility and the effects of predators and productivity on the outcome of trophic interactions

Am. Nat., 1989 134 922 949

http://dx.doi.org/10.1086/285022

Leibold

M.A.

A graphical model of keystone predators in food webs: trophic regulation of abundance, incidence, and diversity patterns in communities

Am. Nat. 1996 147 784 812

http://dx.doi.org/10.1086/285879

Leising

Horner

Pierson

The balance between microzooplankton grazing and phytoplankton growth in a highly productive estuarine fjord

Progr. Oceanogr., 2005 67 366 383

http://dx.doi.org/10.1016/j.pocean.2005.09.007

Litchman

Klausmeier

Trait-based community ecology of phytoplankton

Annu. Rev. Ecol. Evol. Syst., 2008 39 615 639

http://dx.doi.org/10.1146/annurev.ecolsys.39.110707.173549

Maestrini

Granéli

Environmental conditions and ecophysiological mechanisms which led to the 1988 Chrysochromulina polylepis bloom: an hypothesis

Ocean. Acta, 1991 14 397 413

Margalef

Estrada

Blasco

Functional morphology of organisms involved in red tides, as adapted do cecaying turbulence

Taylor

D.L.

Sliger

H.H.

Toxic dinoagellate blooms

North Holland

Elsevier

315 320

May

R.M.

Stability and Complexity in Model Ecosystems 1974

Princenton University Press

Mitra

Flynn

Predator-prey interactions: is 'ecological stoichiometry' sufficient when good food goes bad?

J. Plankton Res., 2005 27 393 399

http://dx.doi.org/10.1093/plankt/fbi022

Mitra

Flynn

Promotion of harmful algal blooms by zooplankton predatory activity

Biol. Lett., 2006 2 194 197

http://dx.doi.org/10.1098/rsbl.2006.0447

Murray

Mathematical Biology 1989

Berlin

Springer

http://dx.doi.org/10.1007/978-3-662-08539-4

Pätsch

Radach

. Long-term simulation of the eutrophication of the North Sea: temporal development of nutrients, chlorophyll and primary production in comparison to observations

J. Sea Res. 1997 38 275 310

http://dx.doi.org/10.1016/S1385-1101(97)00051-8

Pitchford

Brindley

Iron limitation, grazing pressure and oceanic high nutrient-low chlorophyll (hnlc) regions

J. Plankton Res. 1999 21 525 547

http://dx.doi.org/10.1093/plankt/21.3.525

Radach

Pätsch

Climatological annual cycles of nutrients and chlorophyll in the North Sea

J. Sea Res. 1997 38 231 248

http://dx.doi.org/10.1016/S1385-1101(97)00048-8

Sailley

Vogt

Doney

Comparing food web structures and dynamics across a suite of global marine ecosystem models

Ecol. Modell. 2013 261-262 43 57

http://dx.doi.org/10.1016/j.ecolmodel.2013.04.006

Schaffer

W.M.

Ecological Abstraction: The Consequences of Reduced Dimensionality in Ecological Models

Ecol. Monogr., 1981 51 383 401

Smayda

T.J.

Reynolds

C.S.

Community assembly in marine phytoplankton: application of recent models to harmful dinoagellate blooms

J. Plankton Res., 2001 23 447 461

http://dx.doi.org/10.1093/plankt/23.5.447

Solé

Estrada

García-Ladona

Biological controls of Harmful Algal Blooms: A modelling study

J. Mar. Syst. 2006 61 165 179

http://dx.doi.org/10.1016/j.jmarsys.2005.06.00

Solé

García-Ladona

Estrada

The role of selective predation in Harmful Algal Blooms

J. Mar. Syst. 2006 62 46 64

http://dx.doi.org/10.1016/j.jmarsys.2006.04.002

Steele

The structure of marine ecosystems 1974

Harvard University Press

http://dx.doi.org/10.4159/harvard.9780674592513

Teramoto

Kawasaki

Shigesada

Switching effect of predation on competitive prey species

J. Theor. Biol. 1979 79 303 315

http://dx.doi.org/10.1016/0022-5193(79)90348-5

Tillmann

Hesse

Colijn

Planktonic primary production in the German Wadden Sea

J. Plankton Res., 2000 22 1253 1276

http://dx.doi.org/10.1093/plankt/22.7.1253

Touzet

Franco

Raine

Influence of inorganic nutrition on growth and PSP toxin production of Alexandrium minutum (Dinophyceae) from Cork Harbour, Ireland

Toxicon 2007 50 106 119

http://dx.doi.org/10.1016/j.toxicon.2007.03.001

Truscott

Enviromental forcing of simple plankton models

J. Plankton Res. 1995 17 2207 2232

http://dx.doi.org/10.1093/plankt/17.12.2207

Truscott

Brindley

Ocean plankton populations as excitable media

Bull. Math. Biol. 1994 56 981 998

http://dx.doi.org/10.1007/BF02458277

van Donk

Defenses in phytoplankton against grazing induced by nutrient limitation, uv-b stress and infochemicals

Aq. Ecol. 1997 31 53 58

http://dx.doi.org/10.1023/A:1009951622185

Vichi

Pinardi

Masina

A generalized model of pelagic biogeochemistry for the global ocean ecosystem. part I: Theory

J. Mar. Syst. 2007 64 89 109

http://dx.doi.org/10.1016/j.jmarsys.2006.03.006

Wyatt

Horwood

Model wich generates red tides

Nature 1973 244 238 240

http://dx.doi.org/10.1038/244238a0

Yoshida

Hairston

N.H.

Ellner

Evolutionary trade-off between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris

Proc. R. Soc. London 2004 271 1947 1953

http://dx.doi.org/10.1098/rspb.2004.2818