Use of Agent-Based Modeling To Explore the Mechanisms of Intracellular Phosphorus Heterogeneity in Cultured Phytoplankton

There can be significant intraspecific individual-level heterogeneity in the intracellular P of phytoplankton, which can affect the population-level growth rate. Several mechanisms can create this heterogeneity, including phenotypic variability in various physiological functions (e.g., nutrient uptake rate). Here, we use modeling to explore the contribution of various mechanisms to the heterogeneity in phytoplankton grown in a laboratory culture. An agent-based model simulates individual cells and their intracellular P. Heterogeneity is introduced by randomizing parameters (e.g., maximum uptake rate) of daughter cells at division. The model was calibrated to observations of the P quota of individual cells of the centric diatom Thalassiosira pseudonana, which were obtained using synchrotron X-ray fluorescence (SXRF). A number of simulations, with individual mechanisms of heterogeneity turned off, then were performed. Comparison of the coefficient of variation (CV) of these and the baseline simulation (i.e., all mechanisms turned on) provides an estimate of the relative contribution of these mechanisms. The results show that the mechanism with the largest contribution to variability is the parameter characterizing the maximum intracellular P, which, when removed, results in a CV of 0.21 compared to a CV of 0.37 with all mechanisms turned on. This suggests that nutrient/element storage capabilities/mechanisms are important determinants of intrapopulation heterogeneity.

P hytoplankton play an important role in the ecology and biogeochemistry of freshwater and marine systems. Thus, understanding their composition and physiology (e.g., carbon fixation rate) is important. The elemental composition of phytoplankton was originally considered to be relatively uniform (1). However, subsequent research has shown that it can vary between species (2) or clones (3) and within a species under different conditions (4,5). The composition can also vary within a species under the same conditions. Two measurements of intracellular P are commonly made: cell-based P (mol P/cell), which we refer to as P quota, and biomass-normalized P (mol P/mol C), which we call internal P content. While P quota is generally easier to measure on a cellspecific basis, it is the internal P content that is important for many physiological processes. Recent observations of intracellular P in individual cells have shown significant intraspecific heterogeneity in laboratory and field populations, with coefficients of variations (CV) ranging from 0.08 to 2.11 for P quota and 0.10 to 2.10 for internal P content (Table 1).
Heterogeneity in intracellular P may have important consequences because it can result in heterogeneity in physiology (e.g., photosynthesis rate), which can have significant effects on community ecology (6). Specifically, for phytoplankton the specific growth rate (biomass accumulation) is generally considered to depend on the cell quota of the limiting nutrient according to Droop's function (4). Because the function is nonlinear, heterogeneity in cell quota can affect the population growth rate. This is a case of Jensen's inequality, which states that for a concave function, the expected value of the function is less than the function of the expected value (7). If the population's average growth rate is determined based on the population's average cell quota (i.e., average the quota then determine the growth rate for the population), then it will be higher than if the growth rate is determined for each individual and then averaged to get the population-averaged growth rate (8,9).
There are several mechanisms that can create heterogeneity in the intracellular P of a population of phytoplankton. Microscale patches of extracellular (bioavailable) nutrients are created by zooplankton excretion or other processes, and cells randomly encounter these patches, which leads to heterogeneity (10,11). Macroscale mixing (advection and dispersion) causes populations at one location to be made up of cells with different life histories (i.e., exposure to external nutrient concentration), resulting in different levels of intracellular P (9). Genetic differences, when populations are composed of different species and strains or clones, can lead to phenotypic heterogeneity (12). Cell cycle asynchronicity, when populations consist of cells at different stages in the cell cycle and of different sizes, translates to heterogeneity in the P quota (mol P/cell). The internal P content (mol P/mol C) is not directly affected by this mechanism, but it can be indirectly affected if P and C uptake are asynchronous. The cell cycle can also affect nutrient uptake (13). Phenotypic variability, when cell division creates two daughter cells that are not exactly the same in properties (e.g., size) (14) or physiology (e.g., nutrient uptake) (15), leads to heterogeneity in intracellular P directly (i.e., uneven division of P) or indirectly (e.g., as a result of different uptake rates). Feedback between all of these mechanisms can modify the heterogeneity. Intracellular P can affect nutrient uptake, excretion, and mortality (16)(17)(18)(19). Intracellular P also affects the growth rate and cell cycle (4,20).
In the laboratory conditions explored here, asynchronicity in cell cycle, phenotypic variability, and feedbacks between these mechanisms can create heterogeneity. Other mechanisms do not apply to the laboratory. Microscale patchiness is not important, because there are no zooplankton. Macroscale mixing is assumed not to apply to the mixed reactor. Genetic differences apply to populations with multiple strains, which is not the case in a laboratory single-strain culture.
Several studies have shown there can be significant intraspecific heterogeneity in intracellular P (Table 1); theory suggests that it is important and that there are multiple mechanisms that can cause it. However, the importance or relative contribution of these mechanisms and the effect the resulting heterogeneity has on population growth rates are not well understood.
Agent-based modeling (ABM), also known as individualbased modeling (IBM), can be used to simulate population heterogeneity and investigate the mechanisms underlying it. A previous study used this approach to explore the heterogeneity in a field population (10). That study found that microscale patchiness was the main cause of heterogeneity under the field conditions tested. Although Bucci et al. (10) simulated laboratory conditions, they did not compare these results to observations. There is a need to validate model-predicted intracellular P heterogeneity against observation at the laboratory scale. Here, we use the same general approach as Bucci et al. (10) to investigate the heterogeneity observed in laboratory culture. We develop a model which includes a number of mechanisms of heterogeneity and apply it to the data of Nuñez-Milland et al. (21), who grew the marine diatom Thalassiosira pseudonana in batch culture and measured P quota using synchrotron X-ray fluorescence (SXRF). We calibrated the model to the observations. We then performed simulations with different mechanisms turned off to estimate their contribution to the heterogeneity.

MATERIALS AND METHODS
Overview. The ABM simulates one type of phytoplankton (Thalassiosira pseudonana) in culture. Individual cell state variables include the cell size (m; in mol C/cell) and internal P content (q; in mol P/mol C) ( Table 2).
The P quota (mol P/cell) is the product of the cell size and internal P content. Individual cells in the model can grow, take up and excrete nutrient, divide, and die. The growth rate is modeled as a function of photosynthesis and respiration. Photosynthesis is assumed to depend on the internal P content and light, with no effect from temperature and other nutrients. Model cells take up and excrete P. External P is assumed to be in the form of inorganic phosphorus (PO 4 ) (i.e., organic P is not explicitly considered). The model includes mortality, and dead cells are tracked. Details of the processes, modeling approach, and state variables and calculated values used in the model are provided below. The model code is provided in the supplemental material.
Photosynthesis, respiration, and growth rate. Between divisions, a cell's size changes due to photosynthesis and respiration (10): where m is the cell mass (mol C/cell), P (1/day) is the photosynthesis rate, R (1/day) is the respiration rate, and G (1/day) is the growth  rate. The photosynthesis rate is a function of light and the internal P content (4,9,22): where P,MAX (1/day) is the maximum photosynthesis rate, L I is the light limitation term (light [1] or dark [0]), q 0 (mol P/mol C) is the subsistence internal P content, and q (mol P/mol C) is the internal P content. The original formulation proposed by Droop (4) uses the P quota (q and q 0 in units of P/cell), but it is also commonly applied with internal P content (q and q 0 in units of mol P/mol C) (9,23). If q drops below q 0 , then the cell does not perform photosynthesis. Cell division. Cell division is based on the size state variable (m) and a specified minimum cell size parameter (m 0 ). When a cell grows to twice the minimum, it divides and creates two daughter cells with approximately half the size (9,22). Diatoms have a silicate outer cell wall, called a frustule, which provides protection from predators (24). In many diatoms, the frustule makes the cell division process complex, because the daughter cell's valve forms inside the mother cell; therefore, it must have a smaller diameter (25). This gradually reduces the average diameter of the population over time (26). The diameter is increased periodically through the formation of an auxospore (a phase of sexual reproduction) that does not have a rigid silicate shell (26,27). This size reduction mechanism has been modeled before for Cyclotella meneghiniana (10). However, size reduction has not been observed for T. pseudonana. Paasche (28) grew T. pseudonana in culture, and there were no significant changes in the mean cell diameter. Hildebrand et al. (29) also did not observe size reduction in this species and provided evidence that this is because the girdle-band portion of the cell is expandable (allows the area where the valves are created to be wider than the ends of the cells). The cell size can also change due to other factors, like nutrient status (30), but our model simulates nutrient-replete batch culture. Therefore, the model includes the increase due to growth and a drop at division, but the gradual size reduction mechanism or other factors affecting the cell size are not considered.
Phenotypic variability. The model includes phenotypic variability by randomizing the parameters of daughter cells at division. For example, the maximum uptake rate (V A,MAX ) can be varied to account for heterogeneity in nutrient uptake (15). The mass split fraction (Sf-m) and P split fraction (Sf-q-m) can be varied from the ideal value of 0.5 to account for unequal division (14). Here, 12 parameters are randomized ( Table 3). The diameter (d) and density (ñ c ) were omitted, because variability in these parameters would lead to cell surface area heterogeneity, which is already introduced by the mass split fraction. The randomization is done using the function proposed by Kreft et al. (10,31). A value is drawn from a modified normal distribution truncated to lower and upper bounds. The function limits the value to within two standard deviations of the means. Additional upper and lower limits can be specified to avoid unrealistic values (e.g., a split fraction that is less than zero or greater than one). The purpose of the model is to explore nonheritable phenotypic heterogeneity, so a global means is used rather than the mother's value (8). The same general approach, but with a distribution based on the mother's value, could be used to simulate evolution. This would require using a mutation frequency (i.e., it would not necessarily occur at every division event). Further, the parameters would have to be properly constrained to account for physiological trade-offs (32) and avoid the unrealistic evolution of a "Darwinian demon," which can maximize all aspects of fitness simultaneously (33). For example, a mother cell with a higher P,MAX would grow faster, its daughter cells will have a higher average P,MAX and will grow faster, and so on, until the population is dominated by cells with unreasonably high growth rates. Evolution is not expected to be significant during the batch culture experiments investigated here. Considering a reasonable mutation rate (34) and making very conservative assumptions that all mutations result in a phenotypic difference that will affect P quota, less than 3% of the population will have differences in P quota due to this modeled evolution. This cannot explain the observed heterogeneity.

Mortality.
Mortality of phytoplankton has been observed to be a function of environmental stress (18). This was shown in T. pseudonana by Franklin et al. (12), where up to 25% of the cells in a population died after 3 weeks of nutrient limitation. Death is simulated stochastically based on a mortality rate, which increases as the cells' internal P content decreases: where D (1/day) is the death rate, D,BASE (1/day) is the base death rate, and q D (mol P/mol C) is the internal P content at death ( D ϭ D,BASE when q ϭ q D ). The model simulates mortality as a discrete event where each cell has a probability of dying at each time step ( D ⌬t). A random number (between 0 and 1) is drawn and compared to the probability. If the random number is smaller than the probability, the cell dies. Dead cells excrete nutrient but do not perform uptake, photosynthesis, or respiration. Internal P content, uptake, and excretion. An individual cell's internal P content changes due to uptake, excretion, and growth dilution (10) according to the following equation: where V (mol P/mol C/day) is the uptake rate, W (mol P/mol C/day) is the excretion rate, and G q accounts for growth dilution. The uptake formulation is based on Michaelis-Menten kinetics, applied on a surface area basis with an additional term from Thingstad (35) to limit uptake at high internal P contents: where V A,MAX (mol P/m 2 /day) is the maximum specific uptake rate per area, A (m 2 ) is the surface area, K M (mol P/m 3 ) is the half-saturation constant, PO 4 (mol P/m 3 ) is the extracellular concentration of phosphate, and q MAX (mol P/mol C) is the maximum internal P content. The reactor is assumed to be completely mixed, and the concentration of extracellular nutrient is assumed to be uniform. The q MAX parameter specifies the maximum possible internal P content a cell can have under optimal (e.g., high PO 4 ) conditions. As q approaches q MAX , the uptake rate reduces to zero. The surface area is calculated from the cell size (m), assuming a cylindrical shape and specified cell diameter (d) and density (ñ c ). Excretion is a first-order function of the internal P content (10,36): where k w (1/day) is the excretion constant. Excretion is assumed to be in the form of phosphate, the same form as extracellular nutrient. Extracellular phosphate. The concentration of extracellular phosphate in the reactor changes due to uptake and excretion by the phytoplankton cells (10): where Vol (m 3 ) is the volume of the reactor. S R is the superindividual number (see the section below on agent accounting). Batch culture and dilution. The model simulates a batch culture. To model a cell population through growth and subsequent reinoculation (see "Initial and boundary conditions"), a dilution step is included. The dilution adds extracellular phosphate and removes cells. The removal of cells can be simulated in two ways (8). One approach is stochastic, where each agent has a probability of being removed from the dilution, which is equal to the dilution ratio. The alternative approach is to decrease the S R (see the section below on agent accounting) of the agents by the dilution ratio but keep all of the agents in the model. Here, dilution is modeled using the latter approach.
Agent accounting. The very large number of cells in the culture (up to about 2.5 billion in the application presented here) makes it impossible to explicitly simulate each one. For this case, using a 3.8-GHz processor, it would take 32 years to perform one simulation with that many agents. It would also require 100 terabytes of RAM. Therefore, the model uses superindividuals, where one agent in the model represents many individual cells in the real culture (10,37). Each agent has an S R , which specifies the number of real cells it represents. For example, a population of 10 7 real cells can be simulated with a population of 10 3 superindividuals, each representing 10 4 real cells (S R ). The model has two groups of superindividuals: living and dead cells. Each group has a range for its number of superindividuals. The simulation starts with a number of superindividuals at the midpoint of that range. During the simulation, to keep the number of agents within these specified ranges, the model splits or combines cells. If the number of agents goes below the minimum, the model splits the cell with the largest S R into two cells, each with an S R that is half that of the original cell. The two cells retain the other state variables of the original cell. If the number of agents exceeds the maximum, then the model combines the two cells with the lowest S R . The cell's combined S R is the sum of the individual cells, and its size and internal P content is determined by a weighted average of the two original cells to conserve mass. The number of superindividuals used in the simulation is important. When it is too low, binning cells into superindividuals leads to the same problem encountered with the population-level approach (Jensen's inequality). The number has to be sufficiently high to resolve the heterogeneity of the population (i.e., it is a sampling problem). This is established by increasing the number of superindividuals until the results are no longer affected by further increases.  70 Perry (30) Implementation. The model is implemented in the Netlogo ABM framework, version 5.0 (38). A constant time step is used. For each time step, first every agent calculates its uptake and excretion and updates its internal P content; this is done without changing the concentration of extracellular P. Second, after all agents have finished updating their internal P content, the model calculates the new concentration of extracellular P. The agents then will calculate their growth rate, divide, if applicable, and run the mortality routine. If any of the earlier steps causes the number of agents to be outside the acceptable ranges, then the split/combine routine is run. PLM. The model includes an equivalent population-level model (PLM) for comparison (39). The PLM is for the same reactor, with two groups of cells: alive and dead. The PLM does not consider heterogeneity within these groups, which means that each cell within the alive or dead group has the same mass and internal P content and, for the living cells, the same growth rate.
Experimental data. The model is applied to the laboratory experiments of Nuñez-Milland et al. (21), who grew T. pseudonana in batch culture and measured the P quota of individual cells using SXRF. Three cultures were run for different methods of measuring the P quota of cells. Here, the data from the cultures corresponding to the nickel and gold grids are used. Thalassiosira pseudonana (CCMP 1335) was grown in batch culture in a 1-liter vessel. Cells were grown in medium with 1.5 M phosphate, which is a nutrient-replete condition. The cultures were grown with a 14-h light/10-h dark cycle and were stirred manually daily. Cells were sampled during the exponential growth phase (around 11 h into the light phase on the third day) and analyzed using SXRF. The precision of SXRF analysis (quantified as CV) is better than 5% for most elements (40), which is much less than the cell-to-cell variability (Table 1 and Results).
Simulations performed. A number of simulations were performed. The first one, the baseline simulation (M1), is calibrated to the observations. This simulation has all of the mechanisms of heterogeneity turned on, and the variability for each parameter was optimized by minimizing the error between model and data (see the section on optimization results below).
Several simulations then were performed to investigate the contribu-tion of various mechanisms of heterogeneity to the observed heterogeneity (Table 4). Simulations X1 to X12 are the same as the baseline simulation but with one mechanism of heterogeneity switched off. For example, the X1 run has variability in P,MAX turned off. The contribution of each mechanism is determined by comparing the CV of the X run to that of the M1 run. This approach is similar in spirit, but different in two details, to that used by Bucci et al. (10). First, they grouped parameters (i.e., randomization for several parameters was turned on/off together in one simulation), while here individual parameters are considered. Second, Bucci et al. (10) investigated the contribution of mechanisms by performing simulations with only one mechanism turned on. Here, the contribution of mechanisms is investigated using simulations where individual mechanisms are turned off. This was done because some of the mechanisms do not, by themselves, create enough heterogeneity to desynchronize the population. The use of runs that remove a single mechanism allows for determining the contribution of a single mechanism while keeping the population unsynchronized. Initial and boundary conditions. The model starts with a uniform population of living cells (i.e., same m and q). If a mechanism of heterogeneity is turned on, then the cells start with that parameter randomized. The concentration of extracellular P in the media (at the start and for each dilution) is 1.5 M. The model is allowed to spin up for 4 days, after which the population is completely desynchronized in all simulations. A 1:100 dilution then is made (i.e., 99% cells are removed and medium is replaced with new medium). This point corresponds to the beginning of the experiment, and the cell density matches the one measured in the laboratory. The simulation is then run for 3.5 days, stopping 11 h into the light cycle. This is consistent with the experimental procedure where the samples were taken during the exponential growth phase. The run time is sufficient to remove the effects of the initial conditions and to reach a stable P quota distribution (about 8 generations). The ranges of superindividuals for alive and dead cells are 1,000 to 2,000 and 100 to 200, respectively. The time step used is 0.0005 days (43 s).
Parameter values: means and variabilities. For each phytoplankton parameter, means and variability (CV) need to be specified (see the section above on phenotypic variability). The means for the parameters were calibrated manually to match the SXRF data, and most of the resulting  (42). A score of 0 means the error in the X run is the same as the error in baseline, while a negative score means the X run has a higher level of error than the baseline. f The variables in this column are those for which assigned heterogeneity was removed. For example, the X1 run has variability in P,MAX turned off. g This is the observed heterogeneity of the laboratory culture from Nuñez-Milland et al. (21).
values are within the range of the literature ( Table 3). The mean value for the maximum photosynthesis rate ( P,MAX ) is lower than the literature range. This may be because this model does not consider temperature or self-shading limitation. An optimization algorithm was used to find the best variability for each of the 12 parameters (quantified as CV) ( Table 3). The objective of the algorithm is to minimize the model error, which is quantified as the root mean squared error (RMSE) between observed and modeled individual P quotas. The algorithm proceeds in a stepwise manner. At each step, the model is run 10 times (to average out variability between runs) and the RMSE is averaged. At the start, an initial minimum RMSE is determined from the starting conditions. The algorithm then proceeds down the list (randomized) of parameters and increases or decreases (random) the CV by a small amount. The resulting RMSE is compared to the minimum RMSE. If the new RMSE is lower, the change is kept and the minimum RMSE is updated. The algorithm continues until there are no further improvements in RMSE. The parameters for the PLM are set to the means of the corresponding parameters in the ABM.

RESULTS
Optimization results. The optimization algorithm was run for 9,600 model runs (80 passes through 12 parameters with 10 runs each) 9 times (3 different starting points and 3 different random number seeds). The results show that 8 of the runs find a similar final RMSE, while one run gets stuck at a significantly higher RMSE (see Fig. S1 in the supplemental material). This can happen because the model is stochastic and the RMSE will vary among repeat runs. The optimization uses a simple hill-climbing method, so when the algorithm is in a plateau region of the parameter space, the stochastic results can cause the algorithm to get stuck in a phantom local minimum. More sophisticated optimization routines are available (41) which can avoid this problem. However, for the purpose of this investigation, the current routine's results are acceptable. The optimization algorithm found variability for the parameters that was within the range of values used in other models (10,23,31,39) and observed in laboratory and field populations (14,15). The CV values of the optimization run with the lowest error were adopted for the baseline simulation (Table 3).
Model-data comparison. The results from the baseline simulation (M1) are compared to the observations (SXRF) in Fig. 1A, which shows individual P quotas in femtomole P per cell. The data are combined from two replicate cultures (see Materials and Methods). The model reproduces the observations relatively well, with an RMSE of 0.28 (Table 4). There is a group of cells with similar observed P quotas, between 4 and 4.5 fmol P/cell, and this feature is not reproduced by the model. The time from culture inoculation to the SXRF observation was relatively short (3 days), and these cells may represent a synchronized subpopulation. The simulation had a CV for P quota of 0.37, which is close to, but slightly higher than, the data's CV of 0.35.
The model predicts that dead cells account for about 3% of the population, and they are mostly at the lower end of the P distribution (Fig. 1B). The fraction of dead cells is within observed ranges for diatoms in culture. Brussaard et al. (18) found the percentage of dead cells ranges from 2 to 28% depending on the amount of nutrient limitation. The higher percentage of dead cells at the lower end of the P distribution is also consistent with previous observations (12,18). In the model, this is due to a higher death rate for cells with lower internal P content and excretion after cell death.
Sources of heterogeneity. To investigate the contribution of various mechanisms of heterogeneity, a number of simulations with different mechanisms switched off were performed (Fig. 2). The relative contributions were determined by comparing the results of each run with a mechanism removed to that of the baseline simulation. Most simulations have CVs that are similar to the baseline (within the 99% confidence interval) (Fig. 3). However, the simulations that had no variability in q MAX , and to a lesser degree, Sf-q-m produced lower CVs. These two mechanisms produced statistically significant differences in the mean CV value compared to the baseline simulation (Table 4). Most runs had skill   scores close to 0, indicating that their results were similar to the baseline ( Table 4). The run with no variability in q MAX had the lowest skill score, followed by Sf-q-m, which indicates that those runs had substantially more error than the baseline. The skill score (Murphy [42]) compares the mean squared error (MSE) of a reference simulation (baseline) to that of an X run. A score of 0 means the X run has the same skill (error compared to the data) as the baseline, while a negative score means the X run has more error than the baseline. These results show that the main mechanism for creating heterogeneity is variability in maximum internal P content (q MAX ). This parameter specifies the maximum internal P content a cell can have under high external nutrient conditions, which applies to the exponential phase of the batch culture experiments. The variation in the split fraction of intracellular P at cell division (Sf-q-m) makes a smaller contribution. The remaining heterogeneity can be attributed to cell cycle asynchronicity and interactions between different mechanisms. An asynchronous population has heterogeneity in cell size (CV of 0.21), and even if there is no heterogeneity in internal P content (q; mol P/mol C), the P quota (mol P/cell) would be heterogeneous. There are numerous interactions among the various mechanisms. For example, uptake increases q, which in turn increases P , which decreases q (via growth dilution), which then further increases uptake, and so on.
Effect of heterogeneity. Most of the current methods for simulating phytoplankton, as implemented in models with ecological (3,43) or water quality (44,45) objectives, or as described in textbooks (46,47), use a population-level modeling (PLM) approach. In the PLM approach, the population (at one location and time) is assumed to be homogeneous. This assumption can cause a significant error when the growth rate is simulated based on the internal P content using the nonlinear Droop model (9). Specifically, due to Jensen's inequality and the concave shape of the Droop function (equation 2), the population average growth rate, determined by calculating the growth rate of individuals and averaging over the population, is less than the population average growth rate determined by averaging the internal P content of the individuals and then calculating a single growth rate. This is illustrated by comparing the average of P to the P of the average internal P content (see Fig. S2 in the supplemental material). An ABM allows for interpopulation variability and does not incur this error. The effect of this averaging assumption is demonstrated by comparing the results of the ABM to the equivalent PLM. The PLM/ABM growth rate ratio is only 1.03. However, at the end of the simulation, the PLM's population was 40% larger (Fig. 4). The difference between ABM and PLM growth rates is similar to that found by Bucci et al. (10) (ratio of 1.05) for laboratory conditions. The difference in growth rate can be explained by Jensen's inequality, which has a larger effect when heterogeneity is larger. Models of field conditions where additional mechanisms (e.g., microscale patchiness) lead to heterogeneity have higher ratios (1.47 for Bucci et al. [10] and 1.22 for Hellweger and Kianirad [9]).

DISCUSSION
This paper explores the mechanisms that underlie heterogeneity in intracellular P in cultured phytoplankton. An ABM was developed which includes a number of mechanisms of heterogeneity, including phenotypic variability in a number of functions. The model was calibrated for the conditions in the laboratory culture, with parameter averages determined manually and their variation determined by an optimization algorithm. The model was able to reproduce the observed heterogeneity. The laboratory data had a CV for P quota of 0.35, and the baseline simulation (all mechanisms turned on) had a CV of 0.37. Simulations then were performed with individual mechanisms turned off to determine the contributions of individual mechanisms to the overall heterogeneity in the population. The simulations demonstrated that most of the mechanisms have little effect on heterogeneity, only two had a statistically significant difference in median CV value compared to baseline simulation (99% confidence interval), and skill scores were substantially less than 0. The two mechanisms, variation in maximum internal P content (q MAX ) and split fraction of P at division (Sf-q-m), each contributes substantially to heterogeneity, so that turning them off affects the population CV and skill score. Variation in maximum internal P content (q MAX ) is responsible for most of the heterogeneity. Our results suggest that nutrient/ element storage capabilities/mechanisms are important determinants of intrapopulation variability.
The contributions of the various mechanisms are a function of the experimental conditions (e.g., single-strain batch culture un-  der nutrient-replete conditions), and our findings may not be applicable to other scenarios. Contributions of various mechanisms may be different under nutrient-limiting conditions or in natural systems where additional mechanisms contribute to heterogeneity (e.g., microscale patchiness, such as in Bucci et al. [10]). This is supported by smaller variability in heterogeneity under laboratory than field conditions (Table 1). On the other hand, due to the rapid division rates of diatoms, they are often first in successional patterns (48) and may not commonly experience long periods of nutrient limitation in the field (49). This means that the relative contributions determined in this study may also apply to common field conditions.
The conclusion that the variability in the maximum internal P content (q MAX ) is responsible for most of the heterogeneity is biologically plausible. Phytoplankton, including diatoms, can store P in excess of their immediate growth requirements (4,16). This luxury uptake is mostly in the form of polyphosphates. For example, Perry (30) demonstrated that T. pseudonana can accumulate between 19 and 43% of cell P in polyphosphates. Since there is not an absolute requirement for the amount of P taken up by this mechanism, it makes sense that it can be highly variable (thus, q MAX is highly variable). In contrast, for example, there is probably very little a cell could do to change the minimum internal P content. Much of the essential P resides in nucleic acids found in rRNA (3,5), and relatively little variation is possible in this pool. However, recent observations show that phytoplankton, including T. pseudonana, can substitute nonphosphorus (sulfurand nitrogen-containing) lipids for phospholipids under P stress (50). It should be pointed out that if the storage is in the form of polyphosphate, then this should not affect the immediate growth rate. Our model is based on the Droop function, which does not differentiate between forms of P. However, even if the polyphosphate does not affect the immediate growth rate, it will allow for continued growth if the external nutrient is depleted, which could affect population yield and competition between species (51). Storing P as polyphosphates could also affect the variability in Sf-q-m (P split fraction). The polyphosphates are stored as granules (52), and splitting a smaller number of such granules may be subject to stochastic variability.
Heterogeneity in internal P content affects population average growth rate (i.e., Jensen's inequality); therefore, it is important from an ecological perspective. Our results suggest a relatively modest effect (3%). However, the model of Bucci et al. (10) suggests that the laboratory environment is lacking some important mechanisms of heterogeneity (i.e., microscale patchiness), and that the heterogeneity in the field is larger and more important. The literature summary of observations (Table 1) also shows larger heterogeneity in field populations. In addition, recent work on sediment cycling of P has underlined the importance of such polyphosphate pools in biogeochemistry: a significant fraction of P that finds its way into calcium phosphate in sediments is derived from diatom polyphosphates (53), such that variability in these pools could be important to global elemental cycles.
There are several potential model improvements. Comparing model results of cell size to observations and, more generally, another independent data set (validation) would be useful. Also, the current model has only one pool for intracellular P (all of which is available for growth), while in reality intracellular P exists in several different forms (54). The model could be improved by having different pools of intracellular P. For example, the model of John and Flynn (55) tracks polyphosphate and structural and other organic P, with the polyphosphate having no effect on the growth rate.
Future research may build on these results. A more comprehensive picture of heterogeneity in cell composition may be obtained by exploring various fractions of P (e.g., polyphosphates, as described above) and also considering other elements (i.e., carbon, nitrogen, and iron). Single-cell analytical techniques are being applied to understand mechanisms of nutrient storage which control internal cell quotas (51,56), and this information should help to link intracellular heterogeneity to external conditions. Another potential expansion would be to incorporate data on heterogeneity in physiological function. For example, Musat et al. (15) measured inorganic carbon and nitrogen assimilation by individual bacteria. Levy et al. (57) measured the growth rate and gene expression of individual yeast cells. An ideal data set would include measured growth and uptake rates and corresponding internal contents of various forms of nutrients of individual cells. Such data may make possible the development of a more mechanistically detailed model that includes a separate pool for polyphosphate.