## Abstract

We try to answer the question of to what extent details in nutrient uptake and phytoplankton physiology matter for population and community dynamics. To this end, we study how two nutrients interact in limiting phytoplankton growth. A popular formulation uses a product-rule for nutrient uptake, which we compare with that on the basis of synthesizing units. We first fit different nutrient uptake models to a dataset and conclude that the quantitative differences between the models are small. Then we study the sensitivity of phytoplankton growth and zooplankton–phytoplankton interactions (ZPi) models to uptake formulations. Two population models are compared; they are based on different assumptions on the relation between nutrient uptake and phytoplankton growth. We find that the population and community models are sensitive to uptake formulations. According to the uptake formulation used in the ZPi models, qualitative differences can be observed. Indeed, although two models based on functions with similar shapes have close equilibria, these can differ in stability properties. Since stability involves the derivatives of formulas, even if two formulas provide close values, large numerical differences in the stability criterion may occur after derivation. We conclude that mechanistic details can be of importance for community modelling.

## 1. Introduction

In ecosystem modelling, some authors use the simplest possible formulations in their model to describe the minimal number of processes of interest to study the addressed problem. Other authors will use the power of current computers to integrate all the available knowledge concerning the ecosystem, and include as much detail as possible. The effect of details in ecosystem models remains an open issue (see also Nisbet & McCauley 2010), though it has already been discussed in the literature (e.g. Raick *et al.* 2006). The inherent complexity of this question is mainly owing to the fact that, in short, some details play a crucial role at the ecosystem level while others can be omitted. Moreover, the choice of the important details, which is driven by the scientific objectives, should also depend on the ecosystem and the associated environmental constraints.

In this paper, we address the problem of accounting for details in the mathematical formulations of nutrient uptake when dealing with population or community models. We consider this problem in the case of simultaneous limitation by multiple chemical elements. It is now widely accepted that trace metals, such as Fe or Zn, can colimit phytoplankton growth (Hutchins & Bruland 1998; Takeda 1998; Hutchins *et al.* 2001; Quéguiner 2001; Leblanc *et al.* 2005) in addition to macronutrients, such as nitrate and phosphate, and light. We here consider two cases: (i) the case of two limiting macronutrients (e.g. N and P) and (ii) the case where the ability to acquire a macronutrient is reduced when the supply of another nutrient (e.g. a trace metal) is not sufficient. For a long time, the terminology in this area has been confused. We could find either the terms colimitation or multi-limitation for the two aforementioned cases. Now, it seems that the terminology has been clarified, mainly since Saito *et al.* (2008) in which colimitation is clearly defined and subdivided into three categories. According to these authors, the first case we mentioned belongs to ‘independent nutrient colimitation’ (type I colimitation) while the second case is encompassed in ‘biogeochemically dependent colimitation’ (type III colimitation). Note that colimitation type II is associated with ‘biochemical substitution colimitation’ in Saito *et al.* (2008). In the DEB context, we could roughly summarize the previous colimitation denominations as follows: type I colimitation corresponds to the parallel complementary substrate case, type II corresponds to the substitutable substrate case and type III corresponds to the sequential complementary case.

We use the concept of synthesizing units (SU; Kooijman 2010) in order to derive mechanistic uptake formulations. These are defined at the individual level since it corresponds to the appropriate organization level for nutrient uptake. This approach is similar to that proposed in Lorena *et al.* (2010). The formulations for the various types of uptake in terms of SU-dynamics all have the same number of parameters and this number happens to equal that of the usual formulations in the literature. We use a set of data to determine their ability to reproduce measured uptake kinetics. We then include these uptake formulations in population and community models in order to understand how the differences in details at the individual level are transferred to population and community levels. At the individual level, we focus on the uptake mechanisms for the two aforementioned type I and type III colimitation cases.

At the population and community levels, for each choice of uptake formulation, we compare two kinds of population models for phytoplankton, one is an extension of Droop's model (which was developed in the one limiting nutrient situation) while the second is a two reserve–one structure DEB model, based on the SU approach. In fact, when only one reserve and one structure are considered, the Droop model is a special case of the DEB model, in which maintenance is neglected (see Kooijman 2010). In this paper, we consider a two reserves and one structure model for the phytoplankton population. In the extension of the Droop model discussed here, the growth rate is represented by a product of hyperbolic functions of quotas (see §3 for details), since it is a common approach. However, in his own work, Droop (1974) compared a minimum formulation with the standard multiplication of hyperbolic terms and concluded that the minimum formulation was more realistic. Thus in this paper, we refer to the product law growth (PLG) model when we consider the growth model obtained by extending Droop's model. In some sense, the main difference between the PLG model and the DEB growth model comes from the assumption on the link between uptake and growth: the PLG model assumes a direct link between uptake and growth while the DEB model, since it is based on SU mechanisms, contains a more complex relationship between uptake and growth (which is even not explicit, as we will see later). We find that the nature of this relation between uptake and growth can have a strong effect on the results presented in this paper. In the same way, two zooplankton–phytoplankton interaction models are considered (for each uptake formulation), a PLG model and a DEB model.

This paper is organized as follows. In §2, different uptake formulations are compared *per se*, in the situation of colimiting nutrients. In §3, we analyse how the quantitative differences between these formulations, once fitted to a given set of data, may influence phytoplankton population growth dynamics. Then, the effect of the uptake formulation choice on zooplankton–phytoplankton interactions is studied in §4. Our results are discussed in §5 and a conclusion follows.

## 2. Uptake kinetics

In this section, we first recall how the Monod–Michaelis–Menten (M–M–M) formulation can be obtained from an enzyme kinetics mechanism, which is referred to the mechanistic approach in this paper. We extend the approach to multi-substrate situations. Enzyme kinetics are often described either through the concentrations of enzymes or through the concentrations of chemical elements. To deploy a usual approach, we use the concept of concentration to derive the mathematical formulations. However, in DEB theory, the fluxes of elements are preferred. Indeed, the concept of concentration requires a well-mixed environment, which is not encountered in cells. Finally, we introduce the concept of SU, which can be defined as enzymes or complexes of enzymes which do not dissociate from substrates (Kooijman 2010).

### (a) Methodology description in the one limiting nutrient case

When a single nutrient is scarce, the M–M–M formulation is often used to describe the uptake rate of this nutrient. This hyperbolic formula is supported by empirical evidence and has a theoretical basis which can be represented by the following reaction scheme:
2.1
in which *S* is the substrate, *Θ*_{*} the free enzyme catalyst, *Θ*_{S} the enzymatic complex resulting from the association of a free enzyme *Θ*_{*}, a substrate molecule *S* and *P* the reaction product, *b*_{1} the binding rate, and *k*_{1} and *k*_{2} dissociation rates.

The *law of mass action* applied to these enzymatic reactions can be written with the following set of differential equations:
2.2
2.3
2.4
Considering that the total enzyme concentration is constant (*Θ*_{*} + *Θ*_{S} = *Θ*) and that, since enzymatic dynamics are fast, a quasi-steady state is reached (d*Θ*_{S}/d*t* ≈ 0), we obtain:
2.5

This leads to the well-known M–M–M expression given in equation 2.6 (see Murray 1993 or McCarthy 1981 for the full development and corresponding hypothesis):
2.6
where *K* = (*k*_{1} + *k*_{2})/*b*_{1} and *J*_{max} = *Θ**k*_{2}.

### (b) Usual formulations for multi-limiting nutrients

Most ecosystem models including colimitation (generally type I colimitation) currently employ either multiplicative (equation (2.7)) or threshold (equation (2.8)) formulations, which are both empirical formulations:
2.7
and
2.8
where *S*_{j} are the limiting nutrient concentrations and *J*_{max,i} is the maximum uptake rate of nutrient *i*. In this paper, *f*_{j} = *S*_{j}/(*K*_{j,i} + *S*_{j}) where *K*_{j,i} is the half-saturation constant associated with nutrient *j* in the uptake rate of nutrient *i*. Note, however, that *f*_{j} could also be related to intracell characteristics through, for instance, a quota function such as
where *Q*_{j} is called the quota and represents the amount of internal nutrient per cell and *Q*_{min,j} and *Q*_{max,j} are, respectively, the minimal and the maximal amount of quota *Q*_{j} needed for growth (e.g. see Geider *et al.* 1998). In order to be more precise, we recall that in Droop's work, the quota corresponds to chemical elements instead of nutrient, so nitrogen instead of nitrate for instance. If nitrogen in reserve is typically nitrate, we should emphasize that nitrogen in structure is typically built into proteins.

In Saito *et al.* (2008), only these two formulations are discussed in the type I case colimitation based on the studies of Droop (1974) and Rhee (1978) dedicated to the experimental examination of their validity. In both studies, it seemed that the minimum (threshold) law gave the best data fitting but Saito *et al.* (2008) argued that this is not sufficient to conclude, considering that ‘the biochemical response to nutrient stress must result in the up-regulation of high-affinity transport systems even under conditions of multiple nutrient stresses’. A few attempts have been made in proposing mechanistic models for type I (independent) colimitation. Most of them are mentioned in Pahlow & Oschlies (2009). In addition, these authors propose a new phytoplankton optimal growth model (with P, N and light colimitation) which encompasses the observed asymmetrical behaviour (e.g. Healey 1985) of the non-limiting cell quota under N or P limitation. These authors indeed argued that none of the multiplicative and threshold formulations could reproduce this asymmetry since formulations of independent colimitation are symmetrical in essence. This argument does not hold if we consider, for example, that the quota functions *f*_{j} in equation (2.8) may have different expressions for N and P as suggested, for example, in Geider *et al.* (1998). Though interesting, the model of Pahlow & Oschlies (2009) relies on strong hypotheses and is specific to the limiting nutrients N and P. Here, we aim at considering general formulations for type I and type III colimitations *sensu* Saito *et al.* (2008).

As far as we know, few publications deal with models for type III colimitation (e.g. Leynaert *et al.* 2004). Buitenhuis *et al.* (2003) present an empirical model for the zinc–bicarbonate colimitation (Saito *et al.* 2008). Although originally presented as a growth model, nutrient uptake can be derived from it. The uptake flux of nutrient 2 (HCO_{3}^{−} in the original paper) depends on the uptake flux of nutrient 1 (Zn) as
2.9
and
2.10
Equation (2.9) can classically be identified as an M–M–M law in which the saturation constant for nutrient 2, *K̃*_{2}, is a function of nutrient 1 availability and uptake parameters *K*_{1} and *K*_{2}.

### (c) Mechanistic description of colimited uptake

#### (i) General scheme for two limiting nutrients

The mathematical formulations developed in the rest of this section are based on particular cases of the reaction scheme given in figure 1. The differences between type I and III colimitation models are obtained through the assumptions made on the reaction rates (binding and dissociation/release rates).

#### (ii) Type I colimitation

The first case concerns the type I colimitation formulation. We use the scheme provided in figure 1 where we consider the synthesizing units approach. In other words, we assume here that the dissociation rates *k*_{ij} = 0,*i,j* ∈ {1,2}. As in §2*a*, the mass action law can be applied to the reaction scheme of figure 1, which gives
2.11
2.12
2.13
2.14
where
2.15
and *Θ* is the total amount of enzymes, which is constant.

A first simplification of this model can be done by assuming that *b*_{11} = *b*_{12} = *b*_{1} and *b*_{21} = *b*_{22} = *b*_{2} which means that the binding rates for the enzyme's active site and a substrate *i* molecule are the same whether one or the two enzyme-active sites are free. This simplification allows us to reduce the complexity of the resulting mechanistic formulation in order to make it comparable with the empirical formulations.

Furthermore, using the same approach as in §2*a* and given the following definitions: *K*_{1} = *k*/*b*_{1}, *K*_{2} = *k*/*b*_{2} and *J*_{max} = *k**Θ*, the uptake flux of each nutrient is written as
2.16

#### (iii) Type III colimitation

We represent in this section the type III colimitation formulation by considering the scheme in figure 1 and assuming that *k*_{12} = *k*_{21} = 0, *b*_{12} = *b*_{21} = 0. Moreover, we assume that the enzyme is a synthesizing unit, which implies that *k*_{22} = 0. Considering now the resulting reaction scheme, in which the uptake of nutrient 2 is dependent upon nutrient 1 supply (type III colimitation *sensu* Saito *et al.* (2008)), one can write
2.17
2.18
and
2.19

Again, in order to reduce the number of parameters to 3 (as it is in the previous uptake formulations), we assume that *b*_{11} = *A*_{1} *k*_{11} where *A*_{1} = 1 has the dimensions of the concentration of substrate *S*_{1}. In addition, we use the same approach as in §2*a* and given the following definitions: *K*_{i} = *k*/*b*_{ii}, *i* = 1,2 and *J*_{max,2} = *k**Θ*, the uptake flux of nutrient 2 is written as
2.20
2.21
2.22

We also want to highlight here that the formulation given by Buitenhuis *et al.* (2003) can be derived from the same scheme as the formulation we just described above, but with different assumptions. Indeed, let us again consider in the scheme given in figure 1 that *k*_{12} = *k*_{21} = 0, *b*_{12} = *b*_{21} = 0 and *k*_{22} = 0. Moreover, we now assume that *b*_{11}≫*k*, which means that *K*_{1} ≃ 0 with the notation defined in the previous formulations. We then obtain
2.23
2.24
2.25
which indeed corresponds to the same functional form as in the Buitenhuis model (Buitenhuis *et al.* 2003).

#### (iv) Comparison of the different uptake formulations

In this section, we compare five previously mentioned formulations to a set of data. To summarize, they consist of: (i) the product law (PL; equation (2.7)), (ii) the minimum law (min; equation (2.8)), (iii) the mechanistic type I colimitation (MTI; equation (2.16)), (iv) the mechanistic type III colimitation (MTIII; equation (2.20)) and (v) the Buitenhuis model (equation (2.9)). For each formulation, the data fitting is achieved with the least-square method on absolute distance between data and calculated corresponding values. The parameter values obtained for each formulation are provided in table 1⇓⇓⇓.

The dataset concerns enrichment experiments performed on a natural community sampled in the Antarctic ocean, in which diatoms are colimited by iron and silicic acid. The uptake rates of silica for different fixed amounts of iron have been measured. Figure 2 compares the data with the model's output, while figure 3 shows the relative distance between observations and simulations: for each data point, we consider the square of the difference between the data and the calculated corresponding value divided by the data value, and we sum all these ratios. In this example, we see that the mechanistic formulation describing the adequate mechanism (experiments describe a type III colimitation) leads to the best fit. The low performance of the MTI formulation is owing to its low ability to reproduce Si absorption when this nutrient is at low concentration (see figure 2). At low silica concentration (i.e. low *S*_{2}), MT3 becomes linear in *S*_{2}, and responds in a hyperbolic way to *S*_{1}

In contrast, MT1 does not respond to *S*_{1} and is also linear in *S*_{2}
This may explain why colim T1 does not fit the data at low Si, where uptake responds to Fe.

Note that *J*_{max}, the maximum uptake, is the only parameter directly comparable between the five models, and in fact turns out to be the one with similar values in the five models.

However, even if the MT3 formulation is more efficient in this case, the quantitative deviation between the observed data and the five uptake models is rather small in all cases and it follows that the choice of an uptake model based on these data is not obvious. It is then a relevant problem to understand if this choice, when implemented in a population growth model or even in a community model, has an impact on the resulting dynamics. In other words, do the small quantitative differences obtained in the uptake quantification have an impact on phytoplankton growth and on zooplankton–phytoplankton interaction dynamics? Sections 3 and 4 address this question.

## 3. Phytoplankton growth models in a chemostat

### (a) Relationship between uptake and growth

Even in the case where only one nutrient is scarce, the relationship between nutrient uptake and population growth is not simple. The methods used to define this relationship in the literature can be based on a phenomenological approach like in the Monod model, on an empirical approach like in the Droop model or on more or less theoretical considerations like in the Morel or DEB models. In the Monod model, the fundamental assumption is that growth is proportional to nutrient uptake, according to a constant growth efficiency. In the Droop model, a hyperbolic relationship between the internal quota and the population growth rate is assumed, originally supported by empirical evidence in Droop (1958). Morel (1987) provides some empirical relationships between a short-term nutrient uptake and phytoplankton growth on a longer-time scale, in response to variations in external nutrient concentrations. These relationships are based on theoretical support for population dynamics at steady state. However, we argue that in a variable environment, where steady environmental conditions are not likely to occur, there is a need for theoretical considerations specific to a dynamical framework. DEB theory allows such an approach. Here we just consider two models, the Droop model, often used in ecosystem modelling, and a DEB model. Appendix A compares more precisely the models and their properties in the context of DEB theory and in the case of just one limiting resource.

### (b) A Droop growth model extension versus a DEB growth model in a chemostat

In order to compare the phytoplankton growth models, we consider a chemostatic environment. It is a simple situation which, however, allows us to investigate the impact of mathematical formulations on population dynamics. The first model is an extension of the Droop (1958) model given into the situation with two limiting nutrients. The second is the DEB model, which is fully described in Kooijman (2010) for instance. Its structure is illustrated in figure 4. We consider the V1-morph case, which is described by an ODE's system comparable to the extension of Droop's model, called here the PLG model. Indeed, in each model, five state variables are involved. For the PLG model, the state variables are two nutrient concentrations *S*_{1} and *S*_{2}, two intracellular quotas *Q*_{1} and *Q*_{2} associated with the chemical elements of nutrients 1 and 2 and a number of cells *X*. The total phytoplankton biomass in this model is *B*_{PLG} = (*Q*_{1} + *Q*_{2})*X*. The Droop model reads
3.1
3.2
3.3

where *D* is the dilution rate of the chemostat, *S*_{0i} the nutrient concentration in the reservoir of the chemostat, *j*_{Si} the specific assimilation rate of nutrient *S*_{i}, *r*_{max} the maximum specific growth rate of the phytoplankton species and *Q*_{min,i} the minimum quota corresponding to nutrient *i*, where *i* = 1,2.

For the DEB model, there are also two nutrient concentrations *S*_{1} and *S*_{2}, two reserve densities *m*_{E1} and *m*_{E2} and one structural biomass *M*_{V}. The total phytoplankton biomass in the DEB model is *B*_{DEB} = (1 + *m*_{E1} + *m*_{E2})*M*_{V}. The DEB model reads
3.4
3.5
3.6
where *j*_{Ei, A} = *y*_{Ei, Si} *j*_{Si} is the specific assimilation rate, *y*_{Ei, Si} the yield coefficient associated with the transformation of *S*_{i} into *E*_{i} (assumed equal to 1 here), *j*_{EiC} = *m*_{Ei} (*k*_{Ei} − *r*) the catabolic flux from reserve *E*_{i}, *j*_{EiR} = *j*_{EiC} − *j*_{Ei}^{M} − *y*_{EiV} · *j*_{VG} the flux of energy rejected by the growth synthesizing unit to the reserve *E*_{i}, *j*_{Ei}^{M} = min(*j*_{EiM}; *j*_{EiC}) the flux of energy from the reserve *E*_{i} allocated to maintenance, *j*_{VG} = *r* + *j*_{V}^{M} the specific gross growth rate and *j*_{V}^{M} = *j*_{V}^{M1} + *j*_{V}^{M2} the actual flux of energy for maintenance paid from structure where
3.7
3.8
3.9
3.10
This formulation has been proposed in Tolla *et al.* (2007) in order to provide a mechanism for the payment of maintenance on structure when reserve is not sufficient. *ρ*_{V} is the preference parameter. When this parameter vanishes, we get the usual switch model describing the switch from reserve to structure for the payment of maintenance. The interested reader can refer to Kooijman (2010) for more details.

Finally, *r* denotes the specific growth rate. Since the growth of structural biomass is processed by a synthesizing unit from the reserve densities, the specific growth rate is provided by the following implicit relation:
3.11
3.12
Note that in the DEB model (3.4), we do not explicitly specify the fate of excreted reserves. If, for instance, nutrient happens to be a nitrate, then reserve is typically either a nitrate or a protein, the mobilized reserve is then nitrate or amino acids and thus excreted reserve can be nitrate, ammonia or (modified) amino acids. For the sake of simplicity, we just assume that excreted reserves correspond to products which are lost for growth during the experiment time.

The comparison between the DEB and PLG models is made through two approaches. The first one comprises comparison of the mathematical expressions of the models. It is known that in the situation where only one nutrient is limiting (one reserve is sufficient in this case), the Droop model (on which PLG is based) is a particular case of the DEB standard model in which maintenance is neglected. In the present paper, we thus consider a DEB model with two reserves and we assume here that the maintenance is null, in order to see how the PLG formulation may be derived from the DEB model in a similar way as in the single reserve situation. Since in the two reserves DEB model the population growth rate is implicitly given, it is difficult (if not impossible) to find it explicitly and to compare it with the PLG expression. Thus, in order to make this comparison easier, we further assume that the turnover rates of reserves (*k*_{Ei}) are the same for both reserves. Data for colimitation by vitamin B_{12} and P suggest that *k*_{Ei} = *k*_{E}, so this assumption is actually quite realistic. The second approach for the comparison of DEB and PLG is done by using two different but quantitatively close uptake formulations and by analysing how these formulations act on population growth according to the choice of the growth model. For the sake of simplicity, we choose two possible uptake formulations from the five uptake rates described in §2: we compare the PL and the mechanistic type III colimitation (MTIII).

### (c) Derivation of a generalized Droop model from the DEB model

In this section, we perform the first analysis; as explained above, we derive an extension of the Droop model by vanishing the maintenance rates *j*_{EiM} = 0, (*i* = 1,2) and assuming that *k*_{E1} = *k*_{E2} = *k*_{E}. This allows us to simplify the implicit equation (3.11), and we obtain
3.13
where *r*_{max} = 3*k*_{E} and *Q*_{min,i} = 3*y*_{EiV}, while the growth rate used in the PLG model (3.1) reads
3.14
As can be seen, if the Droop is a special case of the DEB model in the one reserve–one structure situation, it is not clear what it is in the multi-reserves situation. Indeed, in this case, the DEB growth rate is not explicit while, by construction, the PLG one is explicit and simply linked to quotas, which themselves are directly related to uptake.

### (d) Effects of the choice of uptake formulations on population growth

To each population growth model, we successively associate two uptake kinetics formulations, namely the mechanistic colimitation type III (MTIII) formulation developed in §2 and the PL formulation. The PL is often used in ecosystem models and the colimitation type III provided the best fit of the dataset in §2*c*(iv). The parameter values used for these formulations are those obtained by the data fitting in §2 (see table 1). In practical applications, when uptake data are available and once the choice of a given uptake formulation has been made, the parameters of this uptake formulation are estimated from the corresponding dataset. In our case, the various models for nutrient uptake fitted the data equally well (except for really low concentrations). Figure 5 illustrates the uptake formulation sensitivity on the phytoplankton growth. The asymptotic cell densities for PL and MTIII uptake are almost identical when implemented in the PLG model, but differed by a factor 7 for the DEB model. It means that the PLG model and the DEB growth model do not respond in the same manner to the choice of uptake formulations, even if these formulations lead to a good fit of a given dataset. This result is discussed in §5.

## 4. Phytoplankton–zooplankton interactions in a chemostat

As we did previously for population growth, we now study the effect of uptake formulation on zooplankton–phytoplankton interaction. More precisely, this section aims at understanding (i) how the effects of uptake formulations highlighted in the previous section are robust in more complex community models and (ii) if the addition of a trophic level enhances or reduces the abovementioned effects. The models used here are the same as in §3 in which we added a state variable *Z* for the zooplankton biomass. In order to keep the model as simple as possible, we do not develop here a full DEB model for the predator but use a Monod model (which is a special case of the DEB model for V1-morphs). In some sense, we assume here that the reserve of the predator has a very fast dynamic and is at equilibrium. We use a Monod formulation for the functional response (also called Holling type II functional response; Holling 1959).

The PLG-Zoo model reads
4.1
4.2
4.3
4.4
where the parameters concerning the nutrients, quotas and phytoplankton growth have been defined in the previous section and *a* is the per capita predation rate, *b* the saturation coefficient, proportional to the handling time, and *e* the conversion efficiency, that is, the amount of zooplankton biomass produced per consumed phytoplankton cell.

The DEB-Zoo model reads
4.5
4.6
4.7
4.8
where the parameters have the same definition and meaning as previously and *y*_{PZ} is the yield coefficient corresponding to the conversion of phytoplankton biomass into zooplankton biomass.

Figure 6 shows that the DEB model is more sensitive to the uptake formulation than the PLG model. Furthermore, this sensitivity is not only quantitative but also qualitative. Indeed, as shown in the graphs in figure 6, if we use the MTIII uptake formulation in the DEB model, the cell density and the zooplankton amount reach constant values while, with the PL uptake rate, they oscillate. Moreover, in the present example, in the former case, the zooplankton is excluded while phytoplankton and zooplankton coexist in the latter case. It follows that ecological interpretation and conclusion may differ drastically according to the choice of formulation. As a consequence, even if an uptake rate formulation fits a set of data rather well, as is the case here with the different formulations (see figure 2), it is not sufficient to validate its use in population, community or ecosystem models. In order to validate a formulation established at a given organization level (e.g. individuals) for its use in a model developed for another organization level (e.g. community), it is thus recommended to check its impacts at each organization level by comparing the corresponding dynamics with data obtained at each level. Finally, it is noticeable that even if the phytoplankton dynamics reach the same equilibrium value for both uptake formulations in the PLG model, the zooplankton amount is quantitatively different (factor 2 in the present example).

## 5. Discussion

We showed in §3 that two different uptake formulations which fit well the same dataset, when implemented in the PLG model, provide rather close dynamics. However, when implemented in the DEB model, they lead to rather different amounts of structure. It means that the DEB model seems to be more sensitive to the uptake formulation than the Droop model. Note that an equilibrium was reached for each simulation and thus the sensitivity observed in the DEB case was quantitative. In §4, we further analysed the sensitivity to uptake formulations in zooplankton–phytoplankton interaction models. Again, when the PLG model is used, we observe a tiny quantitative sensitivity while in the DEB model the sensitivity is stronger and even qualitative: an equilibrium is reached with one uptake formulation and a limit cycle takes place with the other. This sensitivity could be explained as follows. The empirical data only related to constant environments, while the chemostat cultures and community models involved transient kinetics. Population models with similar steady-state dynamics can have different transient dynamics and this transient behaviour can dominate population dynamics. Moreover, phytoplankton–nutrient and zooplankton–phytoplankton interactions can guide the system to very low nutrient concentrations, where the differences between the uptake rates are larger. However, this explanation is not correct. Indeed, in the simulations, we could observe that the nutrient concentrations are never low and even remain in the range where the uptake rates are closer. Actually, the quantitative differences in the uptake not only modify the quantitative values of state variables like cell numbers or structure, but they can also have an impact on stability. The stability of an equilibrium, for example, is determined by the eigenvalues of the Jacobian matrix at equilibrium. The Jacobian matrix involves derivatives of the different equations and then the derivatives of the uptake formulations. Now, even if the different uptake formulations are quantitatively close, the way in which their derivatives are involved in the stability conditions can lead to high quantitative differences. It follows that we can expect different qualitative dynamics. Moreover, this sensitivity is enhanced in models where the relation between nutrient uptake and zooplankton–phytoplankton interaction is complex, which is the case in the DEB model.

## 6. Conclusion

In this paper, we first discussed a mechanistic (SU-based) approach, including detailed descriptions, to derive uptake rate mathematical formulations in the multi-limitation case. We showed, using an example, that such an approach can lead to a satisfying fit of uptake data. Moreover, the mathematical complexity (number of parameters and variables) of these mechanistic uptake formulations is exactly the same as the mathematical complexity of more popular formulations, which permits us to easily integrate them in complex ecosystem models.

In a further step, we investigated the potential effect of the choice of a mathematical formulation among various possibilities for a given process for which data are available. Even if the quantitative difference between the different formulations is small, the implementation of the chosen formula in a model describing other processes can lead to very different quantitative and qualitative dynamics. In some sense, this result extends previous ones, such as Fussmann & Blasius (2005). Fussmann & Blasius (2005) explain why in a simple model the sensitivity of the dynamics to the mathematical formulation of a given process can be important. In their case, they are even able to provide a bifurcation diagram and they show that the structure of this diagram is also sensitive to the choice of the formulation, even if the quantitative differences between the considered formulations are very small. However, this sensitivity might be observed only in such simple models. Our study illustrates, on a more complex model, based on more realistic biological assumptions, that sensitivity to mathematical formulations still occurs. Likewise, Kooijman *et al.* (2004) showed that nutritional details in the growth of zooplankton in a simple closed nutrient–phytoplankton–zooplankton system also had a big effect on the bifurcation diagram for the system, when the total amount of nutrient is used as bifurcation parameter. The nutritional detail was that zooplankton only feeds on phytoplankton structure (as is done in this paper) or also requires phytoplankton reserve; the alga–zooplankton conversion is frequently found to be sensitive to algal phosphorous content in practice. While the steady-state density of zooplankton gradually reduces to zero for declining nutrient amounts in the first situation, zooplankton suddenly goes extinct in the second one at a relatively high amount of nutrient.

In our study, we observed the sensitivity to uptake formulation in the DEB model although it was less obvious in the PLG model. The differences in the population dynamics could be explained by the fact that in the PLG model there is an almost direct relation between growth and uptake. This relation, based on empirical considerations, has been introduced by Droop (1958) in his original paper. However, the relation between nutrient uptake and population growth in the DEB model is much more complex and even not explicit.

Finally, we pointed out the importance of uptake physiological details on population growth and on zooplankton–phytoplankton interactions. It is also important to highlight that the validation of a model using a dataset obtained for an organization level (individuals for instance) is not sufficient to ensure that it can be used in models developed for other organization levels (population or community for instance). In other words, the validation of a model involving processes which take place at different organization levels should be based on datasets obtained for each level.

## Acknowledgements

The authors are grateful to Tânia Sousa and Tiago Domingos and three anonymous reviewers for useful comments which helped to vastly improve the manuscript.

## Appendix A. Comparison of growth models in the context of DEB theory

This appendix aims to compare various popular growth models in the DEB context and with one limiting nutrient/chemical element. Even if it is not directly the topic of the article, which mainly deals with colimitation, it can be helpful to understand how these different growth models behave since, in some sense, parts of this paper extend the results presented here. We here use the standard DEB notations, thus *S* is replaced by *X* (see table 4 for notation details).

#### (a) Andersen model

Andersen (1997) proposed the algal growth model
A 1
A 2
A 3
where
and
where the symbols are given in table 4. Nutrients not only enter cells at rate *v*_{i} but also leave cells at rate *v*_{e}; this leak is not further specified, but serves to motivate the existence of the parameter *S*^{′}. Andersen mentions the problem this model gives for low *S* and high *Q*, which is ‘solved’ by replacing *μ* by max(0,*μ*) for growth and neglecting the effect of leaking nutrients on *S* in these situations.

Like Droop, as we will see later on, Andersen considers simple nutrients (such as phosphate and nitrate) and follows in fact chemical elements, neglecting any overheads. In DEB terms this means that *y*_{EX} = 1, *κ* = 1 and *y*_{VE} = 1/*n*_{EV} where *X* and *E* now stand for some chemical element. Table 4 gives the link with DEB theory, where *m*_{Em} and *ṙ*_{m} now only mean the maximum of *m*_{E} and *ṙ*, and the expressions found for them no longer apply. Even if we can interpret that *v*_{e} relates to maintenance losses, in fact, Andersen thought about an analogy with one-compartment kinetics. Problems with negative uptake and growth rates are avoided assuming that *v*_{e} = 0 and *S*^{′} = 0.

The model in DEB notation reads for *j*_{EM} = 0 and *ḟ*_{m} = *Ḟ*_{m} /*M*_{V} and *j*_{EA} = *y*_{EX} *j*_{XA}
and
where *m*_{Em} is a parameter, like *ḟ*_{m}, *k̇*_{E} and *y*_{EV}. We also have *ṙ*_{m}*k̇*_{E} /(1 + *y*_{EV}/*m*_{Em}).

At steady state we must have for *y*_{XE} = 1/*y*_{EX}

The specific nutrient uptake at constant *X* amounts for *j*_{EAm} = *y*_{EX} *j*_{XAm} to
where *j*_{XAm} is the maximum of *j*_{XA}^{*} as a function of *X*. To avoid the use of *m*_{Em}, we can also write
for half-saturation constant *K* = *j*_{XAm}/*ḟ*_{m} and scaled functional response *f* = *X*/(*X* + *K*).

The parameters are *ḟ*_{m}, *j*_{XAm}, *k̇*_{E}, *y*_{EV} and *y*_{EX} = 1. The latter parameter is included for dimensional purposes. The relationship *j*_{XA}^{*} (*S*^{′}) = *j*_{XM} or *j*_{XM}^{−1} = *j*_{XAm}^{−1} + (*ḟ*_{m} *S*^{′})^{−1} gives perhaps the best map of Andersen's parameter *S*^{′} to DEB concepts, but this map is not free of problems (note that *j*_{XM} = *j*_{EM} *y*_{XE}).

In summary, Andersen's model only deviates from the DEB model for V1-morphs in the uptake rate, which depends on reserve in his model, and in the way maintenance is implemented, by subtracting the costs from assimilation and avoiding the inherent problems. In Andersen's model, the reserve, growth and steady-state dynamics are identical to that of DEB theory. If *m*_{E} is replaced by *m*_{E}^{*} in the expression for *j*_{XA} (*X,m*_{E}), Andersen's model is the DEB model, apart from the overheads and maintenance; that is why the deviating uptake term does not affect the steady-state dynamics. The focus on simple nutrients makes *y*_{EX} = 1 (nutrients are internalized, not transformed) and no growth overhead costs are paid and, as a consequence, no nutrients are released to the environment in association with growth. As a consequence Andersen's model inherits most of the nice DEB properties on homeostasis.

#### (b) Morel model

Morel (1987) modelled a fast short-term uptake (at maximum specific rate *sj*_{XAm} for *s* > 1) in combination with a much lower longer-term uptake (at maximum specific rate *j*_{XAm}) in algae by assuming empirically that nutrient uptake decreases linearly with the reserve density such that
for *j*_{XA} = *fj*_{XAm} and *j*_{EA} = *y*_{EX} *j*_{XA} and maximum reserve density *m*_{Em} = *j*_{EAm}/*k̇*_{E}. If nutrients are just internalized, rather than transformed, we typically have *y*_{EX} = 1. For *s* → 1, the standard food intake is recovered. To avoid negative uptake rates we should have
which always holds for *s* > 1. A change in assimilation does not affect the way growth depends on reserve (density), so *ṙ* = (*m*_{E}*k̇*_{E} − *j*_{EM}/*κ*)/(*m*_{E} + *y*_{EV}/*κ*). A problem with Morel's model to include fast short-term uptake is that it modifies the well-tested long-term uptake. A variant of the Morel model that leaves the long-term uptake unaltered is
To avoid negative uptake rates, we must have for *s* > 1

This model allows the excretion of nutrients in response to starvation. The steady-state behaviour of the modified Morel model is identical to that of DEB for V1-morphs and Andersen. The only difference with the DEB model, in transient state, is that the reserve density approaches its equilibrium faster with factor *s*. For *s* → 1, the modified Morel model fully reduces to the DEB model for V1-morphs (also in transient state).

#### (c) Monod model

The Monod (1942) model is a special case of the DEB model, where *k̇*_{E} → ∞ and *j*_{EM} = 0. It is here included for comparative purposes. The maximum growth rate is then *ṙ*_{m} = *y*_{VE} *j*_{EAm} = *y*_{VX} *j*_{XAm}, with *y*_{VX} = *y*_{VE} *y*_{EX} and *y*_{VE} = 1/*y*_{EV} as before. The comparison is perhaps better served if we select the same maximum growth rate for all models and choose for the Monod model *j*^{′}_{XAm} = (*j*_{XAm}^{−1} + *y*_{VX}/k̇_{E})^{−1}. In this way, the Monod model has similar steady-state dynamics as the DEB model with reserve. Note that this also affects the half-saturation constant *K*.

Table 3.1 of Andersen (1997) gives typical values and ranges for phosphate-limited growth of algae *μ*″ = 1.2(0.8,1.8) d^{−1}, *Q*′ = 3.8(2.5,5.2) µg P (mg C)^{−1}, *μ*′/*μ*″ = 1.15(1.11,1.19), *α*′ = 6.5(2.3,18.2) *l* (d mg C)^{−1}. The corresponding DEB parameters are *ḟ*_{m} = 78 (d mM C)^{−1}, *j*_{XAm} = 13 mmol P(d mol C)^{−1}, *k̇*_{E} = 1.38 d^{−1} and *y*_{EV} = 1.5 mmol P(mol C)^{−1}.

Suppose that the nutrient concentration in the environment oscillates sinusoidally with a period of one day, between 0 and *x*_{m} *K*. Simulation studies show that in the neighbourhood of these parameter values, the predicted growth rates hardly differ for *x*_{m} > 2; only the Monod model predicts slower growth. For lower peak values of nutrient concentration, the predicted growth for the Morel model is higher by a factor *s*, approximately; the other three models behave very similarly. The DEB and the Morel model show least relative amplitude in *m*_{E}, the modified Morel model most; the Morel model shows the largest absolute amplitude in *m*_{E}. The relative amplitude of *m*_{E} of Andersen's model hardly depends on the peak nutrient concentration *x*_{m} *K*. Figures 7 and 8 illustrate the result.

The real significance of reserve is demonstrated if maintenance is present and substrate frequently temporarily absent, and in situations of simultaneous nutrient limitation with frequent absence of any of the required nutrients.

## Footnotes

One contribution of 14 to a Theme Issue ‘Developments in dynamic energy budget theory and its applications’.

- © 2010 The Royal Society