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

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).


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 andKlausmeier 2008).Another example would be the metabolic costs associated with the production of morphological defences (van Donk 1997).The importance of 'prey stoichiometry' in predatorprey 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 1 ) and a dominant one (P 2 ), and a zooplankton herbivore (Z) able to graze on both algal species.The general dynamics of this system can be formulated as in Solé et al. (2006b) ( , , ) ( , , ) ( ) (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 1 * , P 2 * , Z * ), in which source and loss rates are balanced so populations do not vary over time: (2) This state satisfies The following relationships can be derived from this system of equations near an equilibrium point: -Relationship I: Taking the partial derivative with respect to P 1 in (2c), we obtain , D 1 and D 2 do not depend on P 1 , implying that grazers control only the growth of P 2 , leaving P 1 unbounded by grazing pressure.However we have assumed that both grazing functions depend on P 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 P dt P dt 0, 0 1 2   A necessary condition for the dominant species, in our case P 2 , to bloom is to have a greater net growth rate than P 1  According to the system (1), we have Taking the derivatives (on P 1 and P 2 ) of this expression and using I we can obtain the relationships: 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: which implies that, near the equilibrium point, when a small perturbation is applied to the system, the grazing function D 1 must grow faster in the P 1 direction (the P 1 axis in the phase space) than the grazing function D 2 in the P 2 direction.
These relationships are general and independent of the particular functions and parameters adopted.In the next section, we will examine whether the system is stable (it remains near equilibrium after a small perturbation) or unstable (it will diverge with time after a perturbation).This aspect is important in our study because blooms are effectively the result of a perturba-tion 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 andBrindley 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 Applying the relationships I we obtain that Z ˆ1 = Z ˆ2 = 0.Then, the characteristic equation of J would be given by |J − λ Id| = 0, where |J − λ Id| means that the determinant of the matrix, Id, is the identity matrix and λ the eigenvalues of the system.Thus we have An equilibrium solution (P 1 * , P 2 * , 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 1 , P 2 , Z)>(0, 0, 0), but remaining near the origin.A bloom will develop when one population explodes and moves away from the origin.
Therefore, it is interesting to know which constraints determine the instability of the system in this region.According to the RH conditions, this will happen when The term involved in the last inequality (18) will always be positive because the three populations are positive and should increase (P ˆ11 >0, P ˆ22 >0, 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 1 and P 2 will start to grow and, as a response, the grazer will tend to increase (also in agreement with the diagnostic of Gentleman et al. (2003) 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 andNeuheimer 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).
Algal growth will be represented by a logistic expression: B i (P i ) = µ i P i − η i P i 2 , where µi is the intrinsic growth rate species i and η i = µ i /κ i is the quotient between µi and the carrying capacity (κ i ) for species i.A linear mortality function E(Z) = δZ will be adopted for zooplankton.Writing it as an equations system, The initial values and parameters chosen here are derived from previous works (Solé et al. 2006a,b).We consider (P 1 0 , P 2 0 , Z 0 )=(0.1 mg C m −3 , 0.1 mg C m −3 , 6 mg C m −3 ) as representative of initial low population densities near the origin.The phytoplankton saturation parameters are equal η 1 =η 2 =0.0001 m 3 (mg C) −1 day −1 and the mortality parameter for zooplankton is δ=0.1 day −1 .The intrinsic growth rate of each species is fixed at µ 1 =0.56 day −1 and µ 2 =0.49day −1 , thus giving the non-dominant character to P 1 with an intrinsic growth rate greater than that of the blooming species P 2 .After this selection of parameters, we derive the values/ranges for m i , a i and k that verify, using Equations 3 and 4 and Inequalities 6, 7 and 8, the relationships I, II, III with the functions shown in Table 1.The resulting expressions are written in Table 2.The calculations were carried out by means of the Matlab © symbolic tools.
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 1 bloom or result in the same growth rate for the two prey populations.

Model simulations
In order to illustrate the population dynamics of the phytoplankton species, we carried out a numerical integration of model ( 19) using two sets of parameters for each grazing function (Table 4), one of them within the range given in Table 3 and triggering the dominance of P 2 over P 1 and the other one outside this range and leading to dominance of P 1 over P 2 .Figure 1 shows the time evolution of the difference P 2 − P 1 .If the param-Table 1. -Selected functional responses for zooplankton feeding, based on Gentleman et al. (2003) and references therein.The parameters mi and ai 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 P1 and P2.The columns of derivatives are as those reported in Gentleman et al. (2003).
-Conditions that parameters of the different grazing functions must satisfy to allow P 2 to initiate a bloom.P 1 0 , P 2 0 are the initial concentrations of P 1 and P 2 .η 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: D i (P 1 , P 2 ) Range of parameter values satisfying Relationships I-III , 20] eters of the grazing function satisfy the relationships I-III, the net growth rate of P 2 is higher than the net growth rate of P 1 (Fig. 1A), while the opposite happens with the other parameter set (Fig. 1B).The two grazing functions show similar qualitative behaviour although in the Sigmoidal case the divergence between the two populations is faster.Due to its higher nonlinearity near zero, the Sigmoidal function amplifies the differences between the two populations more effectively than the Michaelis-Menten.This discrepancy in the modelled phytoplankton population dynamics highlights the importance of the choice of parameters and type of predation function.
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 condi-tion (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).
Function: D i (P 1 , P 2 ) Satisfying Rel.I-III Not satisfying Rel.I-III lower intrinsic growth rate but otherwise with identical characteristics to the autotrophic flagellates defined originally in the model (Flagellates 1 or P 1 ).The herbivore groups include micro-and mesozooplankton and their phytoplankton grazing function is of a 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 µ 1 =2.02 day −1 for P 1 and µ 2 =1.77 day −1 for P 2 .The initial concentrations of P 2 (0.5 mg C m −3 ) were assumed to be much smaller than those of P 1 , (10 mg C m −3 ) and the model was run over two years.The first 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 2 dominance and one with a different set of parameters that fulfilled these conditions during the whole year (Table 4).It must be noted that in these simulations, several parameters, such as the intrinsic growth rates of the microalgae, are functions of temperature (or other variables) and may vary with time.Thus, in principle, the conditions derived in Table 4 should be recalculated for each time step of the model.However, in a first approach, a parameter can be considered as constant for a time interval if its variability during this period is slower than the characteristic time scale of the system (Vichi et al. 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).

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 andBrindley 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 1 , P 2 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.

Fig. 1 .
Fig. 1. -Difference P 2 −P 1 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 1 =P 2 =0.1 mg C m −3 , Z=6 mg C m −3 .A, satisfying Relationships I-III: P 2 exceeds P 1 (initiation of bloom by P2).B, not satisfying Relationships I-III: P 1 exceeds P 2 (initiation of bloom by P 1 ).The same as the previous cases but for initial concentrations: C, P 1 =0 mg C m −3 , P 2 =0.1 mg C m −3 , Z=6 mg C m −3 and D, P 1 =0.1 mg C m −3 , P 2 =0 mg C m −3 , Z=6 mg C m −3 .

Fig. 2 .
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.

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.

Table 4 .
-Functions and parameter values of the grazing functions (see Table