sm80s1-4271

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 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 (, , , ). 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 (, , , ). Prey selection may allow or inhibit the proliferation of certain phytoplankton species, as has been shown for high nutrient-low chlorophyll (HNLC) regions (). 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 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 (). 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 (). 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 (, ). 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 (, ). Another example would be the metabolic costs associated with the production of morphological defences (). The importance of ’prey stoichiometry’ in predator-prey interactions has been well documented (e.g. , , , ). Nutrient stress may also interact with growth rate and toxin production (, , ). 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 , by means of a simple predator-prey model. These authors used a logistic function for algal growth and a Holling II–type grazing response () 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 , and , 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, 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 (,). 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 METHODSTop 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 :

\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)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: 2847-2875..

– 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 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: 981-998.). 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 = (\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 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 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)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: 2847-2875. 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.

RESULTSTop

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)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: 2847-2875. (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: 2847-2875., Gentleman and Neuheimer 2008Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation. J. Plank. Res. 30: 1215-1231.; 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. 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: 2847-2875.).

Table 1. – Selected functional responses for zooplankton feeding, based on 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: 2847-2875. and references therein. The parameters m_i and a_i are characteristic for each resource and k 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 P₁ and P₂. The columns of derivatives are as those reported in 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: 2847-2875..

	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. 2006aSolé J., Estrada M., García-Ladona E. 2006a. Biological controls of Harmful Algal Blooms: A modelling study. J. Mar. Syst. 61: 165-179.,bSolé J., García-Ladona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 46-64.). 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 P₂ to initiate a bloom. P₁⁰, P₂⁰ are the initial concentrations of P₁ and P₂. η_i = µ_i /κ_i is the quotient between the intrinsic growth rate (µ_i) and the carrying capacity (κ_i) corresponding to species i. The parameters m_i and a_i are characteristic of each prey and k 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 µ₁ = 0.56 day⁻¹, µ₂ = 0.49 day⁻¹ η₁ = η₂ = 0.0001 m³ (mg C)⁻¹ day⁻¹ and δ = 0.1 day⁻¹.

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 Table 1) used in the numerical simulations of the three-equation model (19). µ₁=0.56 day⁻¹, µ₂=0.49 day⁻¹; η₁=η₂=0.0001 m³ (mg C)⁻¹ day⁻¹ and δ=0.1 day⁻¹.

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

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

Full size image

Fig. 1. – Difference P₂−P₁ as simulated with the three-equation model using the grazing functions of Table 1 and the parameter values of Table 4. Initial concentrations are P₁=P₂=0.1 mg C m⁻³, Z=6 mg C m⁻³. A, satisfying Relationships I-III: P₂ exceeds P₁ (initiation of bloom by P₂). B, not satisfying Relationships I-III: P₁ exceeds P₂ (initiation of bloom by P₁). The same as the previous cases but for initial concentrations: C, P₁=0 mg C m⁻³, P₂=0.1 mg C m⁻³, Z=6 mg C m⁻³ and D, P₁=0.1 mg C m⁻³, P₂=0 mg C m⁻³, Z=6 mg C m⁻³.

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. 1995Baretta-Bekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363-379., Baretta and Baretta-Bekker 1997Baretta J., Baretta-Bekker J. 1997. Special issue: European regional seas ecosystem model II. J. Sea Res. 38: 169-413., Baretta-Bekker et al. 1995Baretta-Bekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363-379., 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: 231-248., Pätsch and Radach 1997Pätsch J., Radach G. 1997. 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. 38: 275-310.) 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., Baretta-Bekker 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: 173-193.).

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. 2006bSolé J., García-Ladona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 46-64.). 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ía-Ladona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 46-64., 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. 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: 89-109.), 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).

Full size image

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.

DISCUSSIONTop

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 1989Leibold M.A. 1989. Resource edibility and the effects of predators and productivity on the outcome of trophic interactions. Am. Nat. 134: 922-949., Kretzschmar et al. 1993Kretzschmar M., Nisbet R., MacCauley E. 1993. A predator-prey model for zooplankton grazing on competing algal populations. Theor. Popul. Biol. 44: 32-66., 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): 597-608., 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: 2895-2916., 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): 741-771., Grover 1995Grover J.P. 1995. Competition, herbivory, and enrichment: nutrient-based models for edible and inedible plants. Am. Nat. 145(5): 746-774.) and experiments (Bell 2002Bell T. 2002. The ecological consequences of unpalatable prey: phytoplankton response to nutrient and predator additions. Oikos 99(1): 59-68.). 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 1995Truscott J. 1995. Environmental forcing of simple plankton models. J. Plankton Res. 17: 2207-2232., Pitchford and Brindley 1999Pitchford J., Brindley J. 1999. Iron limitation, grazing pressure and oceanic high nutrient-low chlorophyll (HNLC) regions. J. Plankton Res. 21: 525-547.).

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. 1979Teramoto E., Kawasaki K., Shigesada N. 1979. Switching effect of predation on competitive prey species. J. Theor. Biol. 79: 303-315., 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: 2847-2875., Gentleman and Neuheimer 2008Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation. J. Plank. Res. 30: 1215-1231., 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(1-4): 201-213.). 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. 2004Yoshida T., Hairston N.H., Ellner S. 2004. Evolutionary trade-off between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris. Proc. R. Soc. London 271: 1947-1953., 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: 357-372.) and the presence of morphological defences (van Donk 1997van Donk E. 1997. Defenses in phytoplankton against grazing induced by nutrient limitation, UV-B stress and infochemicals. Aq. Ecol. 31: 53-58.). Toxin production by HAB-forming 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: 169-176.) 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 1981Schaffer W.M. 1981. Ecological Abstraction: The Consequences of Reduced Dimensionality in Ecological Models. Ecol. Monogr. 51: 383-401.). 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, 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.

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(1-4): 201-213.
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): 597-608.
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: 2895-2916.
http://dx.doi.org/10.1016/j.dsr2.2003.07.003

Baretta J., Baretta-Bekker J. 1997. Special issue: European regional seas ecosystem model II. J. Sea Res. 38: 169-413.
http://dx.doi.org/10.1016/S1385-1101(97)00054-3

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: 233-246.
http://dx.doi.org/10.1016/0077-7579(95)90047-0

Baretta-Bekker J.G., Baretta J., Rasmussen E. 1995. The microbial food web in the European regional seas ecosystem model. Neth. J. Sea Res. 33: 363-379.
http://dx.doi.org/10.1016/0077-7579(95)90053-5

Bell T. 2002. The ecological consequences of unpalatable prey: phytoplankton response to nutrient and predator additions. Oikos 99(1): 59-68.
http://dx.doi.org/10.1034/j.1600-0706.2002.990106.x

Buskey E., Hyatt C. 1995. Effects of the Texas (USA) ‘brown tide’ alga on planktonic grazers. Mar. Ecol. Prog. Ser. 126: 285-292.
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: 3-13.
http://dx.doi.org/10.1016/j.jmarsys.2012.08.005

Ebenhöh W., Baretta-Bekker 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: 173-193.
http://dx.doi.org/10.1016/S1385-1101(97)00043-9

Flynn K.J. 2002. Toxin production in migrating dinoagellates: a modelling study of PSP producing alexandrium. Harmful Algae, 1: 147-155.
http://dx.doi.org/10.1016/S1568-9883(02)00028-8

Flynn K.J. 2010. Do external resource ratios matter? Implications for modelling eutrophication events and controlling harmful algal blooms. J. Mar. Syst. 83: 170-180.
http://dx.doi.org/10.1016/j.jmarsys.2010.04.007

Flynn K.J., Iringoien X. 2009. Why aldehyde-induced insidious effects cannot be considered as a diatom defence mechanism against copepods. Mar. Ecol. Prog. Ser. 377: 79-89.
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: 357-372.
http://dx.doi.org/10.1016/0022-0981(95)00140-9

Gentleman W.C., Neuheimer A.B. 2008. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation. J. Plank. Res. 30: 1215-1231.
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: 2847-2875.
http://dx.doi.org/10.1016/j.dsr2.2003.07.001

Grover J.P. 1995. Competition, herbivory, and enrichment: nutrient-based models for edible and inedible plants. Am. Nat. 145(5): 746-774.
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: 169-176.
http://dx.doi.org/10.3354/meps225169

Holling C. 1959. Some characteristics of simple types of predation and parasitism. Can. Entom. 91: 385-398.
http://dx.doi.org/10.4039/Ent91385-7

Holt R., Grover J., Tilman D. 1994. Simple rules for interspecific dominance in systems with exploitative and apparent competition. Am. Nat. 144(5): 741-771.
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: 313-321.
http://dx.doi.org/10.1093/plankt/fbi011

Jessup C.M., Bohannan B.J.M. 2008. The shape of an ecological trade-off varies with environment. Ecol. Let. 11: 947-959.
http://dx.doi.org/10.1111/j.1461-0248.2008.01205.x

Kierstead H., Slobodkin L. 1953. The size of water masses containing plankton blooms. J. Mar. Res. 12: 141-147.

Kretzschmar M., Nisbet R., MacCauley E. 1993. A predator-prey model for zooplankton grazing on competing algal populations. Theor. Popul. Biol. 44: 32-66.
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: 922-949.
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: 784-812.
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: 366-383.
http://dx.doi.org/10.1016/j.pocean.2005.09.007

Litchman E., Klausmeier C. 2008. Trait-based community ecology of phytoplankton. Annu. Rev. Ecol. Evol. Syst. 39: 615-639.
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: 397-413.

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. 315-320.

May R.M. 1974. Stability and Complexity in Model Ecosystems. Princeton University Press. 263 pp.

Mitra A., Flynn K. 2005. Predator-prey interactions: is ‘ecological stoichiometry’ sufficient when good food goes bad? J. Plankton Res. 27: 393-399.
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: 194-197.
http://dx.doi.org/10.1098/rsbl.2006.0447

Murray J. 1989. Mathematical Biology. Springer, Berlin.
http://dx.doi.org/10.1007/978-3-662-08539-4

Pätsch J., Radach G. 1997. 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. 38: 275-310.
http://dx.doi.org/10.1016/S1385-1101(97)00051-8

Pitchford J., Brindley J. 1999. Iron limitation, grazing pressure and oceanic high nutrient-low chlorophyll (HNLC) regions. J. Plankton Res. 21: 525-547.
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: 231-248.
http://dx.doi.org/10.1016/S1385-1101(97)00048-8

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. 261-262: 43-57.
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: 383-401.

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): 447-461.
http://dx.doi.org/10.1093/plankt/23.5.447

Solé J., Estrada M., García-Ladona E. 2006a. Biological controls of Harmful Algal Blooms: A modelling study. J. Mar. Syst. 61: 165-179.
http://dx.doi.org/10.1016/j.jmarsys.2005.06.004

Solé J., García-Ladona E., Estrada M. 2006b. The role of selective predation in Harmful Algal Blooms. J. Mar. Syst. 62: 46-64.
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: 303-315.
http://dx.doi.org/10.1016/0022-5193(79)90348-5

Tillmann U., Hesse K., Colijn F. 2000. Planktonic primary production in the German Wadden Sea. J. Plankton Res. 22: 1253-1276.
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: 106-119.
http://dx.doi.org/10.1016/j.toxicon.2007.03.001

Truscott J. 1995. Environmental forcing of simple plankton models. J. Plankton Res. 17: 2207-2232.
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: 981-998.
http://dx.doi.org/10.1007/BF02458277

van Donk E. 1997. Defenses in phytoplankton against grazing induced by nutrient limitation, UV-B stress and infochemicals. Aq. Ecol. 31: 53-58.
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: 89-109.
http://dx.doi.org/10.1016/j.jmarsys.2006.03.006

Wyatt T., Horwood J. 1973. Model which generates red tides. Nature 244: 238-240.
http://dx.doi.org/10.1038/244238a0

Yoshida T., Hairston N.H., Ellner S. 2004. Evolutionary trade-off between defence against grazing and competitive ability in a simple unicellular alga, Chlorella vulgaris. Proc. R. Soc. London 271: 1947-1953.
http://dx.doi.org/10.1098/rspb.2004.2818