A model for variable phytoplankton stoichiometry based on cell protein regulation

The elemental ratios of marine phytoplankton emerge from complex interactions between the biotic and abiotic components of the ocean, and reflect the plastic response of individuals to changes in their environment. The stoichiometry of phytoplankton is, thus, dynamic and dependent on the physiological state of the cell. We present a theoretical model for the dynamics of the carbon, nitrogen and phosphorus contents of a phytoplankton population. By representing the regulatory processes controlling nutrient uptake, and focusing on the relation between nutrient content and protein synthesis, our model qualitatively replicates existing experimental observations for nutrient content and ratios. The population described by our model takes up nutrients in proportions that match the input ratios for a broad range of growth conditions. In addition, there are two zones of single-nutrient limitation separated by a wide zone of colimitation. Within the co-limitation zone, a single point can be identified where nutrients are supplied in an optimal ratio. When different species compete, the existence of a wide co-limitation zone implies a more complex pattern of coexistence and exclusion compared to previous model predictions. However, additional comprehensive laboratory experiments are needed to test our predictions. Our model contributes to the understanding of the global cycles of oceanic nitrogen and phosphorus, as well as the elemental ratios of these nutrients in phytoplankton populations.


Introduction
Marine phytoplankton take up and assimilate inorganic nutrients, thereby altering nutrient ratios in the ocean. Phytoplankton stoichiometry is, in turn, influenced by environmental factors, as individual cells can regulate their element ratios in response to changes in growth conditions (Rhee, 1978;Goldman et al., 1979). This regulatory capability is controlled by the cell's physiological response traits and, therefore, by evolution. Consequently, species are characterized by different stoichiometries (Geider and La Roche, 2002;Klausmeier et al., 2004a), which may contribute to the maintenance of oceanic biodiversity (Göthlich and Oschlies, 2012;Martiny et al., 2013). Thus, the mechanisms underlying phytoplankton ratios of important elements such as carbon, nitrogen and phosphorus are essential for understanding the biogeochemical cycles of these nutrients, and their role in shaping phytoplankton community composition.
Previous laboratory work has focused mainly on the response of cellular contents and stoichiometry to changes in (i) nitrogen and phosphorus inputs and (ii) dilution rates in continuous cultures ( Fig. 1 and Tables 1 and 2). All this experimental work agrees on some important phenomenology; however, the multitude of different species and environmental conditions used in experiments has led to a lack of consensus with respect to some trends (e.g. phosphorus quota vs dilution rate, Fig. 1c) or, especially, the specific shape of the functional dependence. This uncertainty has limited the ability of those experiments to give generalized answers to fundamental questions. For instance, how do the cell contents of the different nutrients (or quotas) interact with each other, and how does that feedback influence growth and stoichiometry regulation? How does stoichiometry regulation influence nutrient limitation, and the competitive abilities of the cell?
Theorists have developed models intended to reproduce experimental results and identify controls on phytoplankton stoichiometry. In accordance with early experiments (Droop, 1968;Rhee, 1978;Terry et al., 1985), one group of models assumes that growth rate, µ, depends exclusively on the most limiting nutrient (Liebig's law) following a hyperbolic function f (Droop's law): µ = min (f (N), f (P)) (Fig. 2a). These "threshold models" impose, thus, a single functional form f for the growth rate that depends exclusively on whichever nutrient limits growth (Legović and Cruzado, 1997;Klausmeier et al., 2004aKlausmeier et al., ,b, 2007. In consequence, models following this approach predict a similar or "symmetric" response of the nitrogen and phosphorus quotas to changes in the dilution rate of simulated chemostats (Fig. 2b). However, a different trend (i.e. asymmetry) for the quotas has been observed in some experiments (Fig. 1b, c) (Elrifi and Turpin, 1985). To reproduce this trend, some threshold models impose interdependency of the nutrient contents by using phenomenological expressions that assume that the acquisition of one nutrient is limited by the other nutrient (Roelke et al., 1999;Bougaran et al., 2010). A second group of models assumes biochemical independence of the nutrients and a common functional form for the response of growth to each nutrient quota; however, growth depends on the product of those functions (e.g. µ = f (N)f (P)), allowing for limitation by multiple nutrients (Terry, 1982;Saito et al., 2008). A last group of models avoids assumptions about the functional dependence of growth rate on quotas by following a more mechanistic approach. In these models, equations account for essential physiological processes through growth, uptake or respiration terms. This group includes simple models devised to reproduce both marine and terrestrial autotroph behavior (Ågren, 2004), as well as more sophisticated ones able to replicate experimental results from different phytoplankton strains (Pahlow and Oschlies, 2009). The latter example, however, assumes optimal nutrient uptake conditions (Smith and Yamanaka, 2007), which implicitly imposes an evolutionary strategy for the phytoplankton species that yields a fitness maximum, regardless of changes in environmental or predatory pressures. Furthermore, this model also uses expressions that impose quota interactions at the level of nutrient uptake and instantaneous cell acclimation.
These efforts to introduce a dependence of nutrient uptake on cell quotas convey that is essential to account for the well-documented influence of protein regulation on cell nutrient content in order to obtain a reliable model description of stoichiometry. For example, phytoplankton show changes in the number of nutrient-uptake proteins according to the environmental nutrient concentration Gotham and Rhee, 1981a,b), which ultimately can be translated into a cell quota dependence for uptake proteins (Morel, 1987). The resource allocation strategy of a cell therefore determines its stoichiometry. Conversely, the response capacity of the cell is determined by the cell phosphorus content, which is required for ribosomes (P-rich molecules) to synthesize uptake proteins. Thus, accounting for this feedback between resource allocation/protein regulation and nutrient content is essential for predicting phytoplankton stoichiometry.
Previous models have included resource allocation in different ways. For instance, (Klausmeier et al., 2004a) distinguish ribosomes involved in protein assembly from proteins and chloroplasts involved in resource acquisition. This model imposes a fixed stoichiometry for each of these types, and allows the cell to allocate resources as constrained by a tradeoff between assembly and acquisition machinery. A similar functional distinction is made in (Pahlow and Oschlies, 2009), although the stoichiometry and trade-offs are less constrained. Models such as these can replicate many of the patterns summarized in Fig. 1 and Tables 1 and 2. However, the use of prescribed phenomenological forms for growth or quota interactions limits their potential to study the feedback between quotas, cell growth, and stoichiometry.
In this paper, we aim to understand how species stoichiometry emerges from the interplay between environment and cell physiology (and, ultimately, evolution). To this end, we build a theoretical model that, in contrast to previous models, incorporates the key dependence between dynamic uptakeprotein regulation and nutrient availability. We avoid using optimality assumptions or imposing explicit dependencies Table 1. Trend (increasing, ↑, or decreasing, ↓, function) and shape reported in the literature for the cell content (quota, Q) of nitrogen and phosphorus when the growth rate, µ, or the dilution rate, w, are varied in laboratory experiments.

Methods
Our model is composed of dynamic equations for the population content of organic carbon, C, nitrogen, N, and phosphorus, P. Positive terms describe factors that increase population levels of these elements, and negative terms account for decreases. We focus our model on regulation of protein production as a key mechanism underlying the dynamics of the population. Thus, we include explicit equations for the regulation of nutrient uptake proteins. This regulation is controlled by two master functions: a function F that represents protein synthesis and the cues that influence it; and a func-tion G that accounts for the availability of ribosomes to carry out protein synthesis. We consider chemostat conditions, which allows the system to reach steady state as cells are washed out during the dilution process. We use these continuous culture conditions to facilitate a better comparison with existing experimental and theoretical work; nonetheless, the model and results can be easily extended to other environmental conditions.
(i) Equations for organic nitrogen and phosphorus: Population N and P increase at rates V N and V P , respectively, due to nutrient uptake, and decrease due to the washout or dilution process of the chemostat (at a rate w) or through other losses, such as leakage, at rates R N and R P , respectively: The population uptake rate for the different nutrients is given by: (see symbols and units in Table A1). B is the number of cells in the population, V max is the maximum uptake rate of a cell, and K is an effective half-saturation constant accounting for a boundary layer in which the local nutrient  Fig. 2. Curves for growth rate µ, quotas Q, and N : P ratio obtained with a threshold model (Klausmeier et al., 2004b). (a) The growth rate depends only on QX in the X-limited regime, and the shape of the function (hyperbolic) does not depend on which element is limiting. (b) The use of an identical functional form for the dependence of growth on the limiting quota imposes a symmetry of the Q versus w plots around the optimal input ratio (the symmetry between N and P quotas mentioned in the text). Curves for growth rate µ, quotas Q, and N : P ratio obtained with a threshold model (Klausmeier et al., 2004b). (a) The growth rate depends only on Q X in the X-limited regime, and the shape of the function (hyperbolic) does not depend on which element is limiting.
The use of an identical functional form for the dependence of growth on the limiting quota imposes a symmetry of the Q versus w plots around the optimal input ratio (the symmetry between N and P quotas mentioned in the text). concentration is smaller than the bulk nutrient concentration (Pasciak and Gavis, 1974;Mierle, 1985) (see derivation in, e.g. Armstrong, 2008;Bonachela et al., 2011): where K N and K P are the standard (i.e. bulk) half-saturation constants. These constants are a composite of fundamental traits such as cell radius and nutrient ion handling time (Aksnes and Egge , 1991;Bonachela et al., 2011) associated with the uptake of N and P, respectively. D N and D P are the diffusivity of nitrogen and phosphorus in the medium, respectively, and r c is the cell radius. Note that the maximum uptake rate is proportional to the total number of uptake proteins that the cell accumulates at its surface (see Appendix), which can change according to environmental conditions. These variations entail changes in K, which accounts for the observed differences in the half-saturation constant of single species subject to variations in environmental nutrient concentrations (Pasciak and Gavis, 1974). In our model, K = K for large nutrient concentrations (in agreement with the definition of K as a bulk variable), and K > K when the nutrient concentration is scarce (Bonachela et al., 2011).
(ii) Equations for the number of uptake proteins: in our model, phytoplankton are able to regulate the number of proteins, n, they allocate for the uptake of the different nutrients. The cell's content of the nutrient (or quota, Q) is the key factor controlling this regulatory process (Morel, 1987;Song and Ward, 2007;Flynn, 2008). As documented in experimental work for the uptake of a single nutrient (either nitrogen or phosphorus), a population shows an increased number of uptake proteins (or population V max ) when the quota is low, compensating for low uptake rates, and a lower V max when the nutrient is abundant Gotham and Rhee, 1981a,b;Riegman and Mur, 1984;Dyhrman and Palenik, 2001). Thus, in oligotrophic conditions the cell allocates more uptake proteins in order to increase its absorbing area and, therefore, the probability of a successful encounter between the scarce nutrient ions and the protein at the cell membrane. On the other hand, the cell may down-regulate the synthesis of these proteins when the internal concentration of the nutrient reaches the storage maximum Q max in order to decrease biosynthesis and maintenance costs. This strategy translates into an effective increase in affinity under oligotrophic conditions, and a J. A. Bonachela et al.: A model for variable phytoplankton stoichiometry 4345 decrease in eutrophic environments (Bonachela et al., 2011). This mechanism can be encoded in a simple way by using the "protein expression" function, F , given by (Bonachela et al., 2011): where k F is a free parameter controlling the shape of the function and Q min is the minimum amount of nutrient required for the cell to grow. We use here a generic sigmoid function; this choice is motivated by the identification of F with processes involving gene expression -traditionally represented by Hill (sigmoid) functions (Alon, 2007). However, the exact mathematical function is not important, as other normalized functions with similar trends do not alter the qualitative behavior of the model (Bonachela et al., 2011). In addition, it is possible to add other (competing) strategies to the regulation, such as up-regulation of protein synthesis when uptake activity is high and down-regulation when it is low. As justified in (Bonachela et al., 2011), these modifications do not alter the qualitative behavior of the number of uptake proteins. In order to constrain uptake protein synthesis when nutrient levels drop below growth minimum requirements, we also impose the condition F = 0 when Q falls below Q min 1 . Likewise, our model represents a dependence of protein expression on phosphorus: as explained above, phosphorus is a major component of ribosomes, essential for the synthesis of any kind of protein, and also required in ATP, ADP, and other forms of energy storage for the cell (see Sterner andElser, 2002, or discussion in Geider andLa Roche, 2002;Bougaran et al., 2010 and references therein). Thus, we include the "protein repression" function, G, which is strictly a function of the internal phosphorus content of the cell, Q P : where k G,1 and k G,2 are shape factors that help establish the range and boundaries of the function. The lower the content of phosphorus, the stronger the influence of the function on the regulation of the synthesis of proteins. Note that the repression function depends only on Q P because we assume that phosphorus (through ribosomal RNA) is the ultimate limiting nutrient for the synthesis of proteins in our model. Any cellular nitrogen will be allocated to proteins as long as Q N is above the minimum required to grow, Q N min .
Finally, we account for surface area constraints on uptake proteins. For both nutrients, the number of uptake proteins, n, depends on the ratio of absorbing area to total area, A rel , according to the Heaviside function, H (see Appendix for mathematical definitions). This constraint introduces an effective competition between the N-uptake and P-uptake proteins for space on the membrane.
Thus, the equations for the number of uptake sites are: where ν X is the maximum number of proteins for the uptake of nutrient X that a cell can synthesize per unit time (see Appendix). Because the change in n with time depends on the quotas, the maximum uptake rate at any time depends on the nutritional history of the cell. Figure 3 shows the dependence of these equations on the quotas of the different nutrients. The use of the expression and repression functions imposes the allocation strategy of the cell under nutrient limitation. When nitrogen is scarce and the cell is N-limited, the cell prioritizes the synthesis of N-uptake proteins in order to increase Q N and, ultimately, increase growth; only when the N-quota falls below the survival threshold of the minimum quota, Q N min , does the cell down-regulate N-protein synthesis. Quotas below Q N min result in decreasing uptake and growth rates and, ultimately, starvation. On the other hand, phosphorus is the key regulator of biosynthesis due to its role in ribosomes and energy reserves; therefore, P scarcity entails a stronger limitation on the synthesis of proteins, even for quotas above the minimum. Conversely, luxury consumption of a nutrient leads to the down-regulation of the synthesis of its uptake proteins, with cell growth being limited only by the availability of the other nutrient. Note that the different dependencies of n N and n P on the quotas break the symmetry expected a priori from Eqs. (1-6) for N and P.
(iii) Equation for Organic Carbon: the dynamics of the population organic C takes into account a photosynthetic term (only source of organic carbon), maintenance costs, and the loss due to dilution: where P max represents the maximum photosynthetic rate, M C is a maintenance rate and R C a respiration rate. The first term accounts for the synthesis of photosynthetic proteins. Following e.g. (Geider et al., 1998), nitrogen quota plays . (b) Flow diagram describing the different physiological processes in the model: Inorganic nutrients ([N] and [P]) are taken up at rates V N and V P , respectively, and assimilated as quotas, Q N and Q P . Quotas influence the synthesis of uptake and photosynthetic proteins through the expression and repression functions, F and G. Photosynthesis facilitates cell growth, which is also affected by maintenance costs. Uptake costs are hereassociated with V N /N and V P /P. Finally, changes in the number of uptake proteins (i.e. in V max ) influence the uptake of the corresponding nutrient.
a main role in the photosynthetic rate because chloroplasts are protein-(i.e. nitrogen-) rich. In our model, the protein expression and repression functions regulate photosynthesis. The F and G functions for photosynthetic proteins (e.g. Ru-BisCO) increase with the internal content of the nutrient in the cell, as large internal levels of N and P are a proxy for favorable growth conditions (see Fig. 4a). The second term accounts for the maintenance cost associated with uptake proteins, which is represented by the nutrient-specific uptake rate 2 . The third term accounts for the respiration associated with other metabolic processes, and the last term is the mortality or dilution term. The use of a constant P max or a linear respiration term is the consequence of simplification, assuming constant and ideal irradiance and temperature conditions for photosynthesis. This representation of C acquisition con-2 Note that both V X and Q X -and, therefore, X itself -are increasing functions of the growth rate, with . (b) Flow diagram describing the different physiological processes in the model: Inorganic nutrients ([N] and [P]) are taken up at rates V N and V P , respectively, and assimilated as quotas, Q N and Q P . Quotas influence the synthesis of uptake and photosynthetic proteins through the expression and repression functions, F and G. Photosynthesis facilitates cell growth, which is also affected by maintenance costs. Uptake costs are hereassociated with V N /N and V P /P. Finally, changes in the number of uptake proteins (i.e. in V max ) influence the uptake of the corresponding nutrient.
ing the different physiological processes in the model: Inorganic nutrients ([N] and [P]) are taken up at rates V N and V P , respectively, and assimilated as quotas, Q N and Q P . Quotas influence the synthesis of uptake and photosynthetic proteins through the expression and repression functions, F and G. Photosynthesis facilitates cell growth, which is also affected by maintenance costs. Uptake costs are here associated with V N /N and V P /P. Finally, changes in the number of uptake proteins (i.e. in V max ) influence the uptake of the corresponding nutrient.
trasts with the explicit terms for nutrient uptake considered in all the equations above. Thus, our model focuses on nutrient acquisition as the main regulator of the population stoichiometry.
The per-capita rate of change for carbon in the population is used to calculate the growth rate, µ = 1 C dC dt + w. Consequently, the organic carbon per cell in the population remains fixed. This allows us to track the number of cells in the population: Biogeosciences, 10, 4341-4356, 2013 www.biogeosciences.net/10/4341/2013/ and, therefore, to switch from a population-dynamics to a biomass-dynamics approach when required. In this paper, we use Eq. (12) only to calculate the growth rate of the population, and focus on B(t) to characterize its abundance. Thus, the behaviors of Q N = N/B and Q P = P/B are qualitatively similar to those of N/C and P/C, respectively.
(iv) Equations for Inorganic Nitrogen and Phosphorus: to monitor changes in the inorganic nutrients present in the chemostat, we account for the inflow of fresh nutrient (the only positive term in the equations), the outflow due to dilution, and the consumption of the nutrient by the population through the uptake term: where [N 0 ] and [P 0 ] are the concentrations of nutrient in the fresh medium that enter the vessel at the dilution rate. A graphical summary of the physiological processes and interactions considered in the model can be found in Fig. 4b.
(v) Simulation Set-Up: we varied the input ratio, [N 0 ] : [P 0 ], and the dilution rate, w, of the chemostat in simulations to study the reaction of phytoplankton to changes in environmental conditions that may alter nutrient availability. We explored the consequences for phytoplankton stoichiometry by integrating numerically the dynamic equations of the population-level model, Eqs. (1-14). Although the model is applicable to any phytoplankton strain, for the sake of concreteness we used data available in the literature to parametrize our population in accordance with a generic Synechococcus species (see actual values in Table A1). For each pair of (w, [N 0 ] : [P 0 ]), we let the system reach stationary state and calculated different observables such as the number of cells, uptake and growth rates, or nutrient content of the population. We modified the input ratio by altering one nutrient input concentration while fixing the other. This procedure allowed us to compare the qualitative behavior of the model with most of the available experimental and theoretical work. Following (Rhee, 1978), we fixed [P 0 ] and varied [N 0 ] in our analysis.

Results
We first explore the behavior of the per-capita nutrient content (or quota, Q) of nitrogen and phosphorus when the relative input concentration of the nutrients is varied. For low relative input nitrogen (i.e. low [N 0 ] : [P 0 ]) Q N remains fixed at a low value that depends on the dilution rate of the chemostat, while phosphorus reaches saturating levels in the cell (Fig. 5a). Cell growth is limited solely by the availability of nitrogen. As nitrogen input increases, population growth drives inorganic phosphorus to lower values. In consequence, phosphorus becomes limiting as well, and the sta-   tionary value of Q N grows with [N 0 ] : [P 0 ] until the realized maximum storage capacity, Q eff N max , is reached. This effective maximum value differs from Q N max due to the presence of the loss rate R N and the down-regulation of n N , which affect growth even for large concentrations of the nutrient (see Eq. 3 under stationary conditions). Q P shows an analogous pattern as the relative phosphorus input increases (i.e. [N 0 ] : [P 0 ]) decreases (Fig. 5a, inset).
Consequently, in the representation of Q N against Q P (Fig. 5b), zones of exclusive limitation by nitrogen are characterized by the vertical sections of the curve, while zones of exclusive limitation by phosphorus constitute the horizontal parts 3 . For the rest of the points on the curve, cells are co-limited by both nutrients. The range of co-limitation decreases as the dilution rate increases.  In this broad zone of co-limitation (ranging from [N 0 : P 0 ] ∼ 20 to [N 0 : P 0 ] ∼ 70 for the smallest w in Fig. 5), the population N : P matches the input ratio (Fig. 6a). Deviations from this relationship occur at low input ratios because only nitrogen is limiting and Q N remains at its effective minimum for a given dilution rate, whereas Q P remains fixed at its effective maximum value (Fig. 5a). The converse is true for low relative phosphorus ratios (high [N 0 ] : [P 0 ]). These deviation patterns are more pronounced as the dilution rate increases. However, there is a single input ratio at which cells are able to incorporate both nutrients in a proportion that exactly matches the input ratio regardless of the value of the dilution rate, resulting in the reduction of both inorganic nutrients to very low concentrations. This is an optimal nutrient ratio, as the cell is able to draw down nutrients to low levels even at high dilution rates. The ratio at which this happens, [N 0 ] : [P 0 ]| opt , can be determined by representing (Q N : Q P ) / ([N 0 ] : [P 0 ]) against the input ratio for different values of the dilution rate; all these curves intersect at the ordinate 1, with the abscissa corresponding to [N 0 ] : [P 0 ]| opt ∼ 26 in this example (Fig. 6b). This optimal nutrient ratio coincides with the cell quota ratio during exponential growth, and can be calculated from its theoretical definition (see Appendix).
Variations in the dilution rate affect nitrogen and phosphorus quotas differently (Fig. 7). Q N is an increasing linear or convex function of the dilution rate, depending on the limiting nutrient; on the other hand, Q P is a more complicated function that shows convexity when phosphorus is limiting and non-monotonicity otherwise, with a range of change wider than that of Q N (Fig. 7a, b). This pattern translates into a growth rate that, for any input ratio, increases as the quota of the limiting nutrient(s) increases, saturating at a maximum value around µ max = µ N max = µ P max = 0.825 ± 0.005 (Fig. 7c, d). The projected quota at zero growth decreases with declining relative supply of each nutrient, until reaching Q eff min . On the other hand, the quota associated with maximum growth increases as the relative input of the corresponding nutrient increases, until it reaches Q eff max . The latter behavior is more marked in the case of phosphorus than for nitrogen.
We also determined the uptake rate of the limiting nutrient, which follows a hyperbolic functional form depending on the external nutrient concentration (Fig. 8). The maximum uptake rate equations, Eqs. (10) and (11), introduce a dependence of V max on environmental conditions (Bonachela et al., 2011). More specifically, the stationary value of V max for the limiting nutrient decreases as the dilution rate increases, reaching a lower plateau at dilution rates close to the maximum growth rate. At their minimum values, V lo max N ∼ 2.2 × 10 −14 mol cell −1 d −1 and V lo max P ∼ 8.5 × 10 −16 mol cell −1 d −1 . Thus, as environmental conditions improve, the cell allocates more resources to growth and less to uptake (see Eqs. 10-12, or Figs. 3 and 4a).

Model validation
If we aspire to understand the interactions between quotas and their effect on cell growth and stoichiometry, our model must be able to predict realistic behavior with as few assumptions as possible. Indeed, our model makes predictions that match qualitatively most of the phenomenology observed experimentally in phytoplankton populations subject to changing input nutrient ratios or dilution rates. The behavior we observe for the quotas, with Q N increasing and Q P decreasing with the input ratio (Fig. 5a) is described in the early work by Rhee (Rhee, 1978) and Terry et al. (Terry et al., 1985), although in those cases maximum storage limits (corresponding to upper parts of the curves in Fig. 5a) were seemingly never reached. In our model, these effective maximum (i.e. saturation) values for Q are reached as a result of the dynamic equation for each nutrient's uptake protein synthesis.  This physiological range for the quotas imposes limits on protein regulation (Eqs. 10-12) and, ultimately, on the element ratios. That is the case with the N : P ratio (Fig. 6a), for which plateaus at both the upper and lower part of the curve are the result of the cell reaching its extreme quota values (Q eff N max and Q eff P min in the case of the upper plateau, and Q eff N min and Q eff P max for the lower one). This is analogous, for instance, to the situation described in Fig. 3 -compilation of algal cultures -and Fig. 4 in (Hall et al., 2005) (see Fig. 1d). On the other hand, our protein regulation mechanism causes an asymmetric reaction of phytoplankton N and P to changes in dilution rate (Fig. 7a, b) through different dependencies of Eqs. (10) and (11) (or dV max N /dt and dV max P /dt, respectively) on Q N and Q P (Fig. 3). Asymmetry, observed repeatedly in experiments (Terry, 1982;Elrifi and Turpin, 1985), is introduced in other models by using terms devised specifically to that end (e.g. Bougaran et al., 2010).
A last important experimental observation is the hyperbolic dependence of growth on the limiting quotas (see Fig. 1 or corresponding references in Table 1). Imposed in threshold (Droop-like) models (Fig. 2a) (Klausmeier et al., 2004a(Klausmeier et al., ,b, 2007Bougaran et al., 2010), this behavior is an emergent property of mechanistic approaches (see Fig. 7c, d here, or (Pahlow and Oschlies, 2009)). Thus, in each zone of single limitation, threshold models show growth curves that depend only on the nutrient quota and not on the input ratio. Moreover, they impose µ N max = µ P max , in accordance with experiments (Rhee, 1978;Elrifi and Turpin, 1985). In our case, (also expected from Pahlow and Oschlies, 2009), there are marked differences between the single-limitation and colimitation zones: in the single-limitation zones we observe   Fig. 7c, d). The equality of µ N max and µ P max results from the dynamics of our model, as do the specific values of µ N max , µ P max , V lo max N , and V lo max P . These values, contrarily to other models, are not imposed a priori in our case; furthermore, they are in reasonable agreement with static maximum uptake rates measured for Synechococcus (Healey, 1985).

Quota interactions
In addition to replicating empirical findings, our model approaches cell quota interactions differently from other models. The interplay between quotas is a key determinant of phytoplankton stoichiometry, as it decides the response of cells to changes in the environment and, ultimately, cellular elemental ratios. Most other models introduce quota interactions through phenomenological expressions relating V max N and V max P to cell quotas. They impose either an independent limitation (see "dynamic uptake version" in Klausmeier et al., 2004b), or different kinds of cross-limitations (Pahlow and Oschlies, 2009;Bougaran et al., 2010). These expressions assume a direct and instantaneous dependence of maximum uptake rates on nutrient quotas. However, there is no known physiological mechanism by which a cell can instantaneously change its uptake potential. Instead, such adjustments occur through allocation processes that alter the synthesis and degradation rates for uptake proteins (Caperon, 1969;Klausmeier et al., 2007;Song and Ward, 2007). We capture these allocation strategies in our model by allowing nutrient quotas to determine the production rate of uptake sites (Eqs. 10-11) rather than the standing pool of uptake sites (which is proportional to V max ). Thus, our model uses the quotas to specify how cellular resources are allocated to N versus P acquisition.
Therefore, nutrient interactions result in our model from the role of each quota in the expression and repression functions. We avoid imposing phenomenological expressions by constructing these regulatory functions from simple biological arguments. Both F and G encode the physiological responses triggered by changes in the quotas of the different nutrients, which determine the rates of resource acquisition for the cell through photosynthesis and nutrient uptake. Equations (10) and (11) adjust the uptake strategy by altering the production of uptake proteins in the cell, which allows for representation of lags in cellular responses. These simple but mechanistic dynamic equations confer plasticity to the population, allowing it to adjust its stoichiometry in a changing environment. Moreover, we avoid using assumptions about the optimal character of uptake or growth (Klausmeier et al., 2004b(Klausmeier et al., , 2007Smith and Yamanaka, 2007;Pahlow and Oschlies, 2009), thus setting phytoplankton cells in a less-constrained evolutionary context at the expense of a larger parameter space.
Models that, on the contrary, assume no interaction between nutrients show a sharp transition between nitrogen and phosphorus limitation, with a single co-limitation point given by [N 0 ] : [P 0 ]| opt (Fig. 2d) (Tilman, 1982). Based on Droop's law, that constant value equals Q N min /Q P min (Rhee and Gotham, 1980;Klausmeier et al., 2007). For cells described by our model, nitrogen and phosphorus are "interacting essential resources" (Tilman, 1982), showing single-limitation zones but also a wide range of co-limitation (Fig. 5). Consequently, scarcity in one nutrient can be partially compensated by an increased concentration of the other nutrient (Tilman, 1982). Moreover, this co-limitation range depends in our model on the environmental conditions (see Fig. 6a). Experimental work has reported a dependence of co-limitation on the environment, particularly the dilution rate (Kunikane et al., 1984;Elrifi and Turpin, 1985;Terry et al., 1985) and irradiance (Healey, 1985). In addition, our model predicts that the optimal ratio, a specific point in the co-limitation region, is [N 0 ] : [P 0 ]| opt = Q eff N max /Q eff P max (see Appendix). The latter expression, shared by the model in (Pahlow and Oschlies, 2009), is not based on any imposed form for the growth rate, and matches the theoretical definition of the optimal ratio introduced in the Results section (ratio shown by the cell during exponential growth). Furthermore, note that the optimal ratio is far from being related to the Redfield ratio (N : P=16) (Klausmeier et al., 2004a), which in our model is just one more point in the co-limitation zone.

Implications of a broad co-limitation region
As explained in Tilman (1982), Rhee (1978), the outcome of competition between species growing on two non-interacting essential resources depends on the relative distance between the respective co-limitation points. Let us consider the case of nitrogen and phosphorus. When two species, A and B (with single transition points (N : P) opt A and (N : P) opt B , respectively), are placed in the same chemostat and these nutrients are considered perfectly essential (like in threshold models), the theoretical outcome can be determined using the distance from each species' stoichiometric ratio to its respective optimum, = [(N : P) opt − (N : P)] /(N : P) opt . If A and B have different signs, each species will be limited by different nutrients and co-existence is possible (see Fig. 9a). If both have the same sign, the species that is closer to its optimum ratio (i.e. with smaller | |) will, in principle, out-compete the other.
As an important consequence of Fig. 5b, we speculate that a wide region of co-limitation will affect the predicted outcome of phytoplankton species competition. Indeed, the potential outcome in our model is more complicated (Fig. 9b). Opportunities for co-existence do still depend on the competing species being limited by different nutrients (i.e. different sign for their ). is still relevant because it indicates which nutrient is more influential to the growth of each species. However, the input ratio at which the quota of a nutrient reaches its maximum (i.e. saturation) is also important, because it determines the zone where the cell is limited exclusively by the other nutrient. Co-limited phytoplankton have competitive advantage over those limited by a single nutrient (Tilman, 1982). Thus, saturation points influence the outcome of competition experiments in our model, introducing competitive exclusion in parts of the diagram where threshold models predict co-existence (see shaded parts in Fig. 9). Furthermore, in our model both consumption ability (or, ultimately, V max ) and saturation points depend on environmen-  In both panels, gray color indicates co-limitation conditions for the non-interacting nutrients case (a) For non-interacting nutrients, each species' optimal ratio [N 0 ] : [P 0 ] opt matches the point where co-limitation occurs; the optimal ratios of the two species divide the phase diagram, delimiting the co-existence (gray) and exclusion (white) zones. (b) In our model, the optimal ratio is one specific point within a wide co-limitation zone, which is delimited by the ratios at which quotas reach their effective maximum values (saturation points). Fig. 9. Different potential outcomes of competition between two generic phytoplankton species, A and B, for non-interacting (a) and interacting (b) nutrients. The labeled zones represent situations in which there is an opportunity for co-existence. In the rest of the areas, the outcome depends on each species' = [(N : P) opt − (N : P)] /(N : P) opt . In both panels, gray color indicates co-limitation conditions for the non-interacting nutrients case (a). For non-interacting nutrients, each species' optimal ratio [N 0 ] : [P 0 ] opt matches the point where co-limitation occurs; the optimal ratios of the two species divide the phase diagram, delimiting the co-existence (gray) and exclusion (white) zones. (b) In our model, the optimal ratio is one specific point within a wide co-limitation zone, which is delimited by the ratios at which quotas reach their effective maximum values (saturation points). tal conditions and, thus, must be accounted for. This is in contrast to descriptions in which nitrogen and phosphorus do not interact, for which only the (fixed) optimal ratios of consumption rates are needed (Tilman, 1982).
Thus, if real phytoplankton are co-limited by nutrients across a range of input ratios, models featuring wide colimitation zones (like ours or the one in Pahlow and Oschlies, 2009) will best describe competitive outcomes. Moreover, the cell plasticity included in our description allows for predictions about species dominance and co-existence in realistic, dynamic environments.

Experimental test
New experiments are needed to confirm the relationship proposed here between quota interactions, growth, uptake and, eventually, stoichiometry, as no clearly unequivocal and comprehensive experimental evidence exists yet. Such an experimental setup would use two phytoplankton strains, first studied independently and later in competition experiments 4 . To ensure the observation of possible wide co-limitation zones, species need to show broad quota ranges, which must be determined along with the maximum values of growth and uptake rates. Next is the determination of the stationary values of the quotas when the input ratio is varied across a wide range of dilution rates including values close to µ max . These experiments should give a clear picture of the different limitation zones, as well as an estimate of the optimal ratio for each species (Figs. 5 and 6). If enough pairs (w, [N 0 ] : [P 0 ]) are explored, the behavior of the populations with changes in the dilution rate can be also reconstructed (Figs. 7 and 8). Complemented with, e.g. measurements of RNA, these experiments would suffice to provide a clear picture of the cellular allocation strategies under different environmental conditions. As a last step, experiments could be performed in which the two populations compete in chemostats at one (or several) selected dilution rates, and for different input nutrient ratios (Fig. 9). The observation of exclusion in an a priori expected coexistence zone would evidence the existence of a broad co-limitation region. On the other hand, optimal ratios matching Q N min /Q P min and framing the co-existence regions would confirm the classical hypothesis (Rhee, 1978).

The role of size
The parameter values used in this paper correspond to a generic Synechococcus phytoplankton strain. This parametrization has allowed us to compare, at least qualitatively, the obtained results with a large variety of experimental articles available in the literature.
A strength of our model is that it is based on physiological mechanisms that could be parametrized for a range of different phytoplankton species, regardless of the size. Cell radius is explicitly present in Eqs. (5) and (6). These expressions are the result of an approximation (Armstrong, 2008;Bonachela et al., 2011) that has proven to work well for both small and large phytoplankton cell sizes (Armstrong, 2008), while the dynamic equations for the number of uptake sites (or V max ) introduce an effective dependence on size intended to compensate for the coarser performance expected for intermediate sizes. For instance, this expression predicts improved uptake abilities for small cells when it is formulated in terms of a normalized maximum uptake rate, V max /r 2 c (see Armstrong (2008) and Fiksen et al. (2013) for a more detailed study on the performance of the "non-regulation" version of these ex-pressions with respect to size). The dynamics for protein regulation facilitate the description of two well-known extreme results using one single expression: if we use Eq. (5) with Eq. (3) (alternatively, Eq. (6) with Eq. (4)), the resulting form for V becomes independent of the number of uptake sites for scarce nutrient conditions 5 , while the large nutrient concentration extreme provides the maximum uptake rate measured in typical (bulk) experiments. These two cases are usually referred to as diffusion and kinetic limitation regimes (see (Bonachela et al., 2011) for a mathematical derivation, and references therein). These limits remain valid regardless of cell size.
In addition, cells of different sizes have different physiological ranges (e.g. different Q max and Q min ) and, therefore, react differently to similar environmental changes. For small cells with high surface area to volume ratios, physiological regulation of uptake protein dynamics should impact nutrient uptake rates through changes in absorbing area (Eqs. (10) and (11)). However, larger organisms tend to become limited by diffusion more easily and, therefore, the number of uptake proteins may become less important (Pasciak and Gavis, 1974;Armstrong, 2008). Nonetheless, uptake regulation has been observed across a wide spectrum of species and sizes. For instance, the focal species in the classic studies by McCarthy and Goldman (McCarthy and Goldman, 1979) or Gotham and Rhee (Gotham and Rhee, 1981a,b) include Ankistrodesmus falcatus, Asterionella formosa, Euglena gracilis, Scenedesmus sp., and Thalassiosira pseudonana; the radii of these species can be up to two orders of magnitude larger than the value we used in this paper (see Table A1).
Uncertainty about the relationship between cell size and uptake protein regulation or cell stoichiometry could be addressed through additional empirical and modeling studies. From the empirical point of view, the experimental test we propose above could be repeated using several species of different sizes. From the theoretical perspective, one possible approach could follow the approach of (Fiksen et al., 2013). This approach would apply two versions of our model -one with protein regulation and one without-across a continuum of cell sizes. Such a study would reveal the dependence of cellular stoichiometry on cell size in both model versions. The comparison of experimental and theoretical studies could provide important insight into protein and stoichiometry regulation across taxa.

Conclusions
Phytoplankton stoichiometry is a dynamic characteristic of the cell. It results from the allocation strategy adopted by the organism under different environmental conditions. Our model of phytoplankton stoichiometry focuses on protein regulation as vehicle for cell acclimation to different environments. As nitrogen and phosphorus are major regulators of protein synthesis, a realistic description of the cell's resource allocation is essential for our model. Considering allocation leads us to introduce the expression and repression functions, which not only allow for the dynamic regulation of proteins but also introduce other essential features such as interactions between quotas and with growth. The dynamic character of these proteins is eventually translated into an effective dependence of V max on nutrient concentration, which allows this model to capture, for instance, the uptake-related phenomenology observed under diffusion versus kinetic limitation for the cell. In addition, these dynamics yield realistic predictions of other important behaviors for the cell.
Our model represents a phytoplankton cell with a stoichiometric phenomenology more complicated than that inferred from the classic Redfield and Droop work, due to the plasticity of the cell in response to changing nutrient availability. However, the cell's physiological ranges (e.g. maximum and minimum values of the quotas) constrain the plasticity of the cell and are ultimately determined by evolutionary processes. Crucially, these physiological ranges scale with cell size, which also induces competition between the N and P uptake proteins for space at the cell membrane. The effects of cell size are particularly relevant for the extent of nutrient co-limitation. Thus, through their physiological ranges, cells are limited in their ability to match the stoichiometry of their environment, especially during blooms (cells growing close to their maximum growth rate) or in cases of nutrient scarcity (vanishing cell growth).
Our dynamic equations for the number of uptake proteins allow our model to predict how cells respond to changes in the environment. Therefore, our model is suited to describe more realistic situations in which competition occurs under changing ocean conditions, for instance due to diurnal or seasonal variation. Our approach may thus help understand the biotic regulation of oceanic nitrogen and phosphorus, as well as the role of stoichiometry in shaping phytoplankton communities. Nonetheless, the experimental test proposed here should be performed in order to further validate our model and significantly advance knowledge of phytoplankton stoichiometry.

Deductions and definitions
The relationship between the number of proteins devoted to the uptake of nutrient X, n X , and the maximum uptake rate for that nutrient, V X max , is easily obtained from the deduction of the Michaelis-Menten functional form for the uptake rate (e.g. (Bonachela et al., 2011)). The deduction, based on the analogy between the uptake process and an enzymatic reac-tion, states that: where k 2 X is the nutrient ion handling rate and N A is Avogadro's number. The same enzymatic analogy establishes that (Bonachela et al., 2011): where K X is the half-saturation constant for nutrient X, D X is its diffusivity in the medium where it is dissolved, and r X is the radius of reactive part of the uptake protein (see Table A1 for units and values). Thus, if ν X is the maximum number of uptake proteins for nutrient X that the cell can synthesize per unit time, the maximum change in V X max per unit time is given by: The handling constant also plays a role in the determination of the total absorbing area for the cell and, therefore, in the ratio absorbing : total area. Assuming, for simplicity, that r N ∼ r P = r s : which is the main argument of the "competition-for-space" term of Eqs. (10) and (11), represented by the Heaviside function: H (x) = 1 for x ≥ 0 (and 0 elsewhere).
Deduction of the optimal ratio: at [N 0 ] : [P 0 ]| opt , the stationary concentrations of both nutrients in the chemostat are negligible. Thus, from Eqs. (13) and (14), V N /V P = [N 0 ]/[P 0 ]. On the other hand, the optimal ratio is the only point where this happens even for the maximum value of the growth rate. At those growth rates, the cell quotas reach their maximum values, and both V N and V P equal their respective maximum uptake rates, which in turn take their lower values, V lo max N and V lo max P . By using the stationary solution to Eqs. (3) and (4), and assuming loss rates much smaller than the dilution rate (or, alternatively, R N ∼ R P ), we find that In our case, Q eff N max /Q eff P max ∼ 26, which is in accordance with the values deduced from Fig. 6b. Table A1. List of variables and parameters used in this manuscript, representing a generic Synechococcus strain -values collected from Healey (1985); Ikeya et al. (1997); Hense and Beckmann (2006); Pahlow and Oschlies (2009);Flynn et al. (2010); Bonachela et al. (2011).

Symbol Description
Units Value