Validation of a Mathematical Model for Green Algae ( Raphidocelis Subcapitata ) Growth and Implications for a Coupled Dynamical System with Daphnia Magna

Toxicity testing in populations probes for responses in demographic variables to anthropogenic or natural chemical changes in the environment. Importantly, these tests are primarily performed on species in isolation of adjacent tropic levels in their ecosystem. The development and validation of coupled species models may aid in predicting adverse outcomes at the ecosystems level. Here, we aim to validate a model for the population dynamics of the green algae Raphidocelis subcapitata, a planktonic species that is often used as a primary food source in toxicity experiments for the fresh water crustacean Daphnia magna. We collected longitudinal data from three replicate population experiments of R. subcapitata. We used this data with statistical model comparison tests and uncertainty quantification techniques to compare the performance of four models: the Logistic model, the Bernoulli model, the Gompertz model, and a discretization of the Logistic model. Overall, our results suggest that the logistic model is the most accurate continuous model for R. subcapitata population growth. We then implement the numerical discretization showing how the continuous logistic model for algae can be coupled to a previously validated discrete-time population model for D. magna.


Introduction
Studies of the population dynamics of phytoplankton and their zooplankton predators in lentic habitats have found a variety of patterns.Plankton communities have been observed to either oscillate in low or high amplitude cycles or to remain relatively stable throughout the year [1].The same lake may exhibit stability on a given year but switch to oscillation during the following year, and vice versa.A variety of explanations have been proposed for this behavior in the field, including predator-prey interactions, temperature fluctuations, and external influences on nutrient content [2][3][4].Fewer studies have attempted to answer the question of what drives these population dynamics in the laboratory setting.Of the predator-prey models that have been proposed for plankton communities, most do not consider certain elements of zooplankton biology such as density-dependent mortality or age-specific fecundity.These traits are crucial for describing the population growth of zooplankton such as Daphnia magna [5] and may be an important factor when modeling population dynamics that are observed in lakes and other ecosystems.
The aim of the present study is to validate and parametrize continuous models for the growth of green algae (Raphidocelis subcapitata) in the absence of predation by zooplankton.The broader goal is to couple a validated green algae model with a validated discrete-time population model for Daphnia magna [5], and ultimately with an analogous continuous-time model, in order to investigate the possibility of oscillations in the laboratory setting similar to those found in lentic environments.In our previous study of D. magna population dynamics [5], we carried out laboratory experiments in which green algae were fed to D. magna populations on a daily basis, and populations were reared in media optimized for daphnia survival, but not necessarily ideal for algae growth.In particular, it was not known whether green algae could proliferate in daphnia media to an extent that would affect daphnia population dynamics.Thus, it is of central importance to quantitate the rate at which green algae grow in daphnia media and whether this growth has the potential to significantly alter the fecundity and survival of daphnia.In theory, such changes would thereby affect the quantification of population level risk assessments involving experimental exposure of daphnia populations to environmental perturbations, e.g., toxins or temperature change.
We tested several commonly used growth models for organisms with a limiting nutrient: the Gompertz, the Logistic (continuous and discretized), and the Bernoulli population models.We note that the Bernoulli model is a generalization of the continuous Logistic model, which allows for nested model comparison.Each model has been used to describe populations across many scenarios associated with saturated growth processes in biology.We collected experimental data from three replicates of green algae grown in isolation of predation.We describe goodness of fit of several mathematical models for green algae growth in the context of asymptotic theory and bootstrapping techniques.We provide estimated parameter values and computed confidence intervals for the model predictions.Finally, we implement a numerical scheme that can be used to approximate the concentration of green algae on a daily time scale in order to combine the continuous green algae model that we had the most confidence in with a discrete time population model for Daphnia magna.We performed simulations of an unvalidated coupled green algae and daphnia model in order to explore the possible effects of green algae growth on daphnia population dynamics.

Data
To observe the growth of Raphidocelis subcapitata populations (previously known as Pseudokirchneriella subcapitata and Selenastrum capricornutum), we seeded three beakers containing 1 L of media reconstituted from deionized water for Daphnia magna culturing (previously described in [6]) and recorded the population density for eight days.Each population was kept in an incubator (Thermo Fisher Scientific, Waltham, MA, USA) at 20 • C on a 16/8-h light/dark cycle (6 AM-10 PM light, 10 PM-6 AM dark).Media lost to evaporation was replaced daily with deionized water in order to retain a 1 L volume and avoid replenishing nutrients.The 1 L algae cultures were uncovered and inspected for contamination during measurements.We selected a seeding concentration of 7 × 10 7 cells based on previous studies of algal growth in order to observe both the early (growth) and late (saturation) stage dynamics of the population [7][8][9].We measured the density of each population replicate twice using a hemocytometer (Hausser Scientific, Horsham, PA, USA) at 9 AM, 3 PM, and 9 PM daily in order to obtain sufficient data points for parameter estimation and uncertainty quantification.The two measurements of density at each time point for each replicate were averaged to minimize human measurement error.This yielded a total of 23 data points for each of the three replicates.

Asymptotic Theory
The goal of this paper is to determine the most accurate model for algae growth in the absence of consumption.The uncertainty in parameters for each model can be quantified via asymptotic theory.In this section we provide the theory behind our asymptotic theory methodology.The estimation of parameters using asymptotic theory requires mathematical models of the form and the corresponding observation process where θ = (q, y 0 ) ∈ R p+ p is the vector of unknown parameters, q is a vector of p model parameters, y 0 is the number p of initial conditions that is unknown, and C maps the model solution y(t, θ) in R l to the observed states f (t, θ).We consider the initial condition to be unknown because of measurement error.In this investigation, the observation operator will always produce a scalar, and thus C maps R l to R. In fact, in all our considerations we have p = l = 1, i.e., the models are scalar and C = I.
Due to the discrete nature of our experimental data, the observations for our statistical error model occur at n = 23 discrete times t i .Thus, the observations will be To account for measurement error, we use the statistical model for our observations, where E i is a zero mean random variable representing identically, independently distributed (i.i.d.) noise that causes our observed data to deviate from our model solution, and θ 0 is the hypothesized "true" or "nominal" parameter vector that generates the observations {Y i } n i=1 .The existence of this "true" parameter vector θ 0 is a standard assumption in frequentist statistical formulations.The i.i.d.nature of the error in our model implies that E(E i ) = 0 for each i, and that E i = 1, . . ., n, are identically distributed with variance σ 2 0 .Since E i is a random variable, Y i is a random variable with corresponding realizations y i .Asymptotic theory seeks to estimate θ 0 by creating a random variable Θ whose realizations for a given data set y i will be estimates θ of θ 0 .These estimates θ will approximate the true parameters θ 0 , and are obtained by minimizing the ordinary least squares (OLS) cost functional [10,11] where Y = (Y 1 , Y 2 , . . ., Y n ) T .Thus, with Ω being the space of admissible parameters and y i being the realizations of the random variable Y i , provides an estimate for θ 0 .The process of estimating parameters from data is known as an inverse problem, and all inverse problems in this experiment are computed using fmincon in MATLAB (Mathworks, 2015b, Natick, MA, USA, 2015) with function and step tolerances of 10 −20 and 1000 iterations.
Once we have an estimate θn OLS , we wish to ascertain the statistical properties of the estimator Θ.Although we do not know the distribution of the estimator Θ n OLS , we can approximate it under asymptotic theory (as n → ∞) by the multivariate Gaussian distribution [10][11][12] Θ n OLS ∼ N(θ 0 , Σ n 0 ) where, based on previous assumptions, the covariance matrix Σ n 0 is approximated by Here χ n is the sensitivity matrix where θ k is the k th component of the vector θ ∈ R 1×p .The unbiased estimator for Note that dy dt = g(t, y(t), θ) is the differential equation for the green algae model and f (t j , θ) = y(t j , θ) is the forward solution of each model.Because we know analytical formulas that provide solutions for dy dt = g(t, y(t), θ), we can solve Equation (11) by setting up a differential equation in terms of the sensitivity [11].
The χ n matrix provides a measure for how sensitive the mathematical model is to each of its parameters.This can be used to estimate the p × p covariance matrix, Σ n 0 , In order to determine the confidence we have in the parameter estimates, we also compute the asymptotic theory based standard error SE( θk ) = Σn kk for the k th parameter.

Boostrapping
We implemented bootstrapping techniques to complement our asymptotic theory approach with regards to estimating parameter uncertainty.We again assume that we are have experimental data for a dynamical system from an underlying observation process in Equation ( 4) where E i are i.i.d. with mean zero and constant variance σ 2 0 and θ 0 is the "true value" hypothesized to exist in statistical treatments of data [10].The random variable also has realizations We use the following algorithm [10] to compute the bootstrapping estimate θBOOT of θ 0 and its empirical distribution.
1. First estimate θ0 from the entire sample {y i } n i=1 using OLS. 2. Using this estimate, define the standardized residuals ri for i = 1, . . ., n, where n is the number of data points, and p are the number of model parameters.Set m = 0, which will represent the total number of artificial samples we will create.3. Create a bootstrapping sample of size n using random sampling with replacement from the data (realizations) {r 1 , . . ., rn } to form a bootstrapping sample {r m 1 , . . ., rm n }.

Create bootstrap sample points
for i = 1, . . ., n. 5. Obtain a new estimate θm+1 from the bootstrapping sample {y m i } using OLS.6. Set m = m + 1 and repeat steps 3-5 until m ≥ 1000 (this can be any large value, but for these experiments we used M = 1000).
We then calculate mean, standard error, and confidence intervals using the formulas: where θ BOOT denotes the bootstrapping estimator.We present the results of these techniques as standard errors about the mean of the parameter estimates, as well as the parameter distributions created.This procedure is performed for each replicate in our experiments.

Model Comparison: Nested Restraint Tests
We used nested model comparison tests to determine the loss in accuracy by constraining certain models, i.e., holding some parameters constant.In general, we assume that we have an inverse problem for the model observations f (t, θ) and are given n observations with the cost function described above in Equation (5).We are interested in using data to question whether the "true" parameter θ 0 can be found in a subset Ω H ⊂ Ω, for which we make the same assumptions as Banks, Hu, and Thompson [10].Thus, we want to test the null hypothesis H 0 :θ 0 ∈ Ω H , or that the constrained model provides an adequate fit to the data.We then define and θn where y is a realization of Y.It is important to note that S n (y; θn H ) ≥ S n (y; θn ).We define the nonnegative test statistics and their realizations, respectively, by and Tn = T n (y) = S n (y; θn H ) − S n (y; θn ).
We refer to [10] for a description of asymptotic convergence as n → ∞, which yields the model comparison result with the corresponding realizations ûn = U n (y) which can be compared to a χ 2 distribution with r degrees of freedom.In this project we use a χ 2 (1) table when comparing the results from the Logistic model to those from the Bernoulli model.

Akaike Information Criterion
In some cases (such as comparison between the Logistic and the Gompertz), the models are not nested (although they are related through a limiting process-see below) and hence we cannot use the model comparison tests outlined above.However, we can use an alternative model evaluation framework and implement the Akaike Information Criterion (AIC c ) with a small size sample correction [10] in the context of an ordinary least squares framework where n is the sample size and p is the number of unknowns (parameters).This will allow us to suggest which model provides a better fit for the data (models with smaller AIC c values provide better fits).While other goodness of fit tests may be useful for selecting models, we chose to use AIC c , since it is a widely adopted measure of model accuracy (see Sections 4.2, 4.3 of [10]).

Logistic Model
The first model we consider is the widely used logistic model for bounded growth of a population P(t), given by the differential equation where R is the intrinsic growth rate, and K is the carrying capacity for the population under consideration.

Bernoulli Model
We also analyze the data within the context of a Bernoulli model due to Richards [13], given by the differential equation The Bernoulli model has three model parameters, R, K, and β.The parameter β extends the logistic model to allow flexibility in the growth dynamics by allowing the inflection point to change while keeping the carrying capacity approximately the same [13].Setting β = 1 yields the logistic model in Equation (26); hence, the logistic model is a special case of the Bernoulli model, thereby enabling us to use nested model comparison tests described above.

Gompertz Growth Model
The next model we consider is the Gompertz growth model, which is widely used for biological and economic phenomena where population growth is not symmetric about the point of inflection, i.e., growth rates are time dependent.The differential equation form of this model is ), where K is the carrying capacity and κ scales the time.For both the Logistic and Gompertz models, we let X 0 represent the initial condition, i.e., P(t 0 ) = X 0 .
The Logistic and Gompertz models, while not nested, are related through a limiting process.Since we find that the Gompertz model is the limit as ν → ∞ of the generalized logistic model for ν > 0 ). (30)

Logistic Model: Numerical Discretization
Another model we consider is a discrete numerical approximation of the Logistic model.We note that the continuous models described above were simulated using the ode45 algorithm in Matlab.In order to ensure that the logistic model can be coupled to a discrete time model for a D. magna population in which the population size is updated once per day [5], we investigated the logistic model using a forward Euler scheme that was discretized on an hour basis.In this paper, we refer to this discrete Euler-method logistic model as the DEL model.This numerically discretized logistic model is given by the difference equation where P t is the population at time t hours and P t+1 is the population at the next time step.The parameters R and K are analogous to those of the continuous Logistic model and can be interpreted as the intrinsic population growth rate and carrying capacity, respectively.We refer to the initial population, P t=0 , as X 0 in our results and data fitting procedure.

Data fitting and Model Comparisons
Overall, we found that all models provide a reasonable fit to the data.Figures 1-4 show results of the least squares estimation for the three different replicates of the Logistic, Bernoulli, DEL, and Gompertz models.These figures contain 68% and 95% confidence bands around the fits to data.These were constructed by generating 1000 random parameter sets from a normal distribution described by the mean and standard error obtained by the asymptotic theory results, computing the model for each of these parameter sets, and then calculating the respective confidence intervals from model generated points f (t i , θ k ), where k = 1, . . ., 1000 [5].One primary difference between the fits to data that we found was that each model tends to underestimate the initial data and the Bernoulli model provided the closest fit (Figure 4).We note that we chose to estimate the initial condition due to measurement error associated with a low cell density as well as how much influence these discrepancies in error would affect the outcome of the model.The results in Table 1 show the small sample corrected Akaike Information Criterion (AIC c ) scores based on Equation (25) for each replicate and each model.These results suggest that the discrete and continuous logistic population models are better able to describe the green algae growth data than either the Gompertz or Bernoulli models, although the differences in AIC c values were small between all models.The results of the model comparison test performed between the continuous Logistic model and the Bernoulli model are given in Table 2.In addition to fixing β to 1 to reduce the Bernoulli model to the Logistic model, we also fixed the initial condition X 0 at the seed value to enunciate the model improvement with an unalterable initial condition.Of all comparisons, only the Logistic model with X 0 fixed on replicates 2 and 3 yielded a significant result at the α = 0.9 confidence level (or χ 2 (1) .9= 2.706).This indicates that, in general, the OLS cost associated with the Bernoulli model was significantly improved by fixing β and reducing it to the Logistic model.However, restricting the Logistic model further by fixing X 0 does not significantly affect the cost of the Bernoulli model.The Logistic model also has benefits with regards to identifiability, which will be seen in subsequent passages of this document.

Uncertainty Analysis
We compared results of parameter estimation and multiple uncertainty quantification techniques (asymptotic theory and bootstrapping) for the Logistic, Gompertz, and Bernoulli growth models, as well as the numerically discretized version of the Logistic model (DEL model), since each of them provided reasonable fits to the data.We first note that the usual assumption of i.i.d.residuals required for uncertainty analysis held for all models investigated (Supplementary Figures S1-S4).Although methods involving autocorrelation on residuals may be used to investigate the i.i.d.assumption, we investigated this assumption by visually inspecting residual plots, since there were not enough data to perform autocorrelation tests.Visual inspection of residual plots is a commonly used procedure when not enough longitudinal data are available (see [10], Section 3.6).We also note that the normality assumption for the parameter distributions in asymptotic theory was confirmed by our bootstrapping results in all but the Bernoulli model (Supplementary Figures S5-S17).We divide our analysis of the results from uncertainty quantification among the sets of parameters with similar interpretations for each model.For example, each model has an initial condition, a growth rate, and a parameter governing the saturation of growth due to population density.

Uncertainty Analysis: Initial Condition
The bootstrapping distribution results of X 0 estimation are presented in Supplementary Figures S5-S8.These appear to be approximately normally distributed, with some exceptions occurring where the estimates are close to the zero boundary.Bootstrapping estimates of uncertainty for X 0 are compared to asymptotic theory in Supplementary Tables S1-S4.We observed that the parameter estimates for X 0 for all replicates were lowest in the Gompertz model and highest in the Bernoulli model.The standard errors also varied between the asymptotic and bootstrapping techniques depending on the model.For example, the order of magnitude for the standard errors from bootstrapping was greater than asymptotic theory for the Gompertz model and the DEL model.

Uncertainty Analysis: Growth Rate
Each model that we investigated has a parameter that describes the population growth rate (Logistic and DEL: R, Gompertz κ, Bernoulli R).Because numerical estimates for the growth rate parameters will not be equal across models, we analyzed their consistency and uncertainty across replicate data sets within the same model.The bootstrapping distributions for the growth rates for each model were normally distributed except for the Bernoulli model, which was skewed to the right (Supplementary Figures S9-S12).We postulate that one reason for this skewness may be that the Bernoulli model is over-parameterized.Similarly, the standard errors computed for the growth rate within each model were of reasonable size and of the same order of magnitude except for the Bernoulli model (Supplementary Tables S5-S8).The growth rate estimates for the logsitic model differed between the continuous version and its numerical discretization using the euler method (DEL model).Specifically, the growth rate estimates for the DEL model were consistently higher.In addition, the asymptotic standard errors for the DEL model were lower than for the continuous logistic model.

Uncertainty Analysis: Saturation Parameter
The saturation parameter K has the same interpretation for all models we investigated, it is the carrying capacity of the green algae population.The estimates for K were remarkably similar across all models (Supplementary Tables S9-S12).The standard errors were inconsistent between asymptotic theory and bootstrapping for the Bernoulli and Gompertz models.For example, the asymptotic standard error for the estimate of K in replicate 1 for the Bernoulli model was 0.1138, whereas the bootstrapping error was 0.0158 (Supplementary Table S12).These results are important, because the asymptotic standard error would result in a much wider confidence band around the model fit to the data, which is indeed the case for replicate 1 of the Bernoulli model (Figure 4, Top).Since the bootstrapping distributions for K for all of the models are normally distributed, this indicates that the bootstrapping standard errors are accurate (Supplementary Figures S13-S16).

Uncertainty Analysis: Bernoulli Model Parameter β
The parameter β is unique to the Bernoulli model, and scales the rate at which the green algae population reaches carrying capacity.In particular, when β = 1, the Bernoulli model reduces to the logistic model.We found that the estimates for β with asymptotic theory and bootstrapping were >1.We can not confidently say that these estimates are accurate, since the corresponding standard errors are unreasonably high for both uncertainty techniques and for all three replicates (Table 3).Moreover, the bootstrapping distributions for β were not normally distributed and heavily skewed to the right for all three replicates (Supplementary Figure S17), indicating the possible presence of correlations with other model parameters.We speculate that the parameters β and K may be correlated for the Bernoulli model.Our uncertainty analysis results indicate that the Logistic model and its numerically discretized counterpart, the DEL model, are the most accurate models among those we investigated.We forgo a summary of the evidence for this conclusion until the discussion section.Our ultimate aim of validating a model for green algae growth was to couple it to a model for D. magna population dynamics.In this section, we couple two validated models for algae and Daphnia population growth to create a theoretical, unvalidated model to explore potentially complex predator-prey dynamics.
The D. magna model we use is one that we recently validated with experimental population data [5].The validated D. magna model is a specification of the following discrete-time discrete-age structured population model: The population is divided into one-day age classes up to some maximum age i max and the number of daphnids of age i at a time t is p(t, i).The average fecundity of each age class i is given by a(t, i) and the survival rate for daphnids of age i is given by b(t, i).The validated functional forms are a(t, i) = α(i)(1 − q) M(t−τ) and b(t, i) is defined piecewise as µ(1 − c) M(t) if i ≤ 4 and µ if i ≥ 5.A summary of the parameters and variables in the model are listed in Table 4 (see [5] for further details).
We coupled the D. magna population model to the DEL green algae model by assuming that D. magna consumes green algae and that the density-dependent survival and fecundity rates of D. magna are influenced directly by the algae concentration.We modeled the algae population with predation as where δ is a predation coefficient, and M t is the total Daphnia biomass at time t.We chose this functional form based on the assumption that Daphnia consume aglae at a rate proportional to the density of aglae and the total biomass of the Daphnia population, as adult daphnids consume food at a higher rate than younger ones.We used a 24 h time discretization to model algae growth for our simulation study and transformed parameters accordingly, setting δ = 0.001, K = 0.4559, R = 1.34, and the initial algae population P 0 = 0.0633.We modeled the algal influence on Daphnia fecundity and survivorship by setting a(t, i) = α(i)(1 − q) 1/P t−τ and b(t, i) = µ(1 − c)  We found that a coupled model could result in both steady state dynamics as well as oscillatory behavior for different choices of parameter values in the Daphnia model (Figure 5).The deterministic simulations in Figure 5 left show steady state behavior with small deviations relative to the population size.The seemingly random perturbations are due to the age-dependency of the density-indpendent fecundity rate α(i).We found that lowering the density-dependent survival competition parameter c yielded sustained oscillations, and increasing it led to both populations reaching a steady state.Other parameter combinations may also yield similar dynamics, but detailing those values is not the aim of the present study.

Discussion
Our results highlight the importance of performing uncertainty quantification in validated biological models, even in the simple case of saturating growth dynamics encountered for green algae.For example, the ordinary least squares regression seemed to indicate that each of the models we investigated provide reasonable fits to the algae growth data.In addition, parameter estimates were consistent between replicate data sets for each model.From this information alone, one might conclude that the Bernoulli model was the best performing model, since it best fit the initial data X 0 (Figure 4).However, a deeper investigation with uncertainty analysis allowed us to generate confidence bands around the fits to the data, showing that the Bernoulli model was the model for which we could have the least confidence.This result was corroborated by a thorough examination of the standard errors for each of the parameters with similar interpretations across all models.For example, the growth rates for the Bernoulli model had unreasonably high standard errors, while the uncertainty in the corresponding growth rates for the other models were relatively low.Our results also indicated that the Gompertz model had inconsistent standard errors between asymptotic theory and bootstrapping for the initial condition and saturation parameter, emphasizing the importance of using multiple uncertainty quantification techniques to ascertain the best validated model.We observed that the logistic and DEL models have different confidence regions for the model solution in Figures 1 and 2. We attribute the change in computed confidence regions to the differences in numerical discretization, the time step for the DEL model was equal to one hour while the logistic model had a coarser mesh.
We collected replicate data and our results had strong agreement across the three replicates.Although methodology exists to fit all three data sets simultaneously to the same model, we chose to fit them separately to test whether the model validation results were consistent across several repeated experiments.We noticed a slight trend in the residuals resulting from fits of the model to data for each replicate: the fit sometimes either over-or under-estimates the data in groups of threes (Supplementary Figures S1-S4).We surmise that this phenomenon may be explained by human measurement error; some people tend to over or under count the algae when using a hemocytometer.Since all models were confounded with this possibility for human error, we can assume that human error did not influence the analysis by favoring one model over any other.In future work, more accurate and consistent cell counts may be performed with a flow cytometer.Alternatively, a spectrophotometer may be used to approximate algae concentrations.
Overall, we suggest that the population growth of Raphidocelis subcapitata is most accurately modeled using the Logistic equation among the simple growth models we investigated.It is important to note that our findings are limited to controlled laboratory settings with unchanging temperature, constant photoperiod, and no change in nutrient availability.For example, we did not consider the possible influence of photoperiod (light/dark) conditions on algae growth.We also did not consider the affect of limiting nutrients such as carbon, nitrogen, phosphorous, or sulfur.The models we investigated here represent our first approximation of algae growth and seemed to fit the data well even without considering how light affects algae growth and that the affect of limiting nutrients could be described by saturating algae growth.The simplifying assumption we made that growth parameters are independent of light conditions may be investigated in future work to yield a closer fit to the data.Our work here serves the purpose of coupling our green algae model with one of zooplankton (D. magna) population growth in a laboratory setting, e.g., for toxicity testing, but should not be directly extrapolated to populations in lakes.In order to develop an accurate model of community fluctuations in the field, we will need to consider predation by various zooplankton and microbes, competition with other algae, nutrient fluctuation, abiotic drivers, and habitat heterogeneity.This study is, however, a useful step toward developing a more comprehensive model.In particular, our results showing that a coupling a validated green algae model with a validated daphnia model are important because it exemplifies the possibility of using a mathematical model to recapitulate the oscillatory dynamics seen in nature.In contrast, the previously validated daphnia model that did not include algae dynamics was not able to produce oscillations and only resulted in steady state behavior [5], a result that did not account for the broader range of plankton dynamics seen in natural systems [1].We note that our ultimate goal is to validate a coupled continuous time daphnia/algae model since continuous time models, such as the Sinko-Streifer model, are described by partial differential equations (PDEs) with a continuously structured variable and can be more amenable to the estimation of age-dependent parameters than a discretely structured model [14][15][16][17].In this work, we investigated the dynamics of a coupled algae/daphnia discrete-time model as a coarse approximation to understand the qualitative impact of including a dynamic food source on daphnia populations.Our finding that increasing the density-dependent survival constant (c) in a coupled predator-prey model yields oscillatory dynamics compliments previous work that has predicted limit cycles based on increased mortality [18].In an ecological setting, changes in the parameter c could reflect differing nutrient (algae) requirements on the density-dependent survival of daphnia due to differences in species size (e.g., Daphnia pulex vs. Daphnia magna) or other increased sources of density-dependent mortality such as predation on daphnia.Changes in c may also be induced toxicologically.For example, endocrine toxins are known to alter the molt cycle of adult daphnids through incomplete ecdysis [19], which may have an indirect affect on density-dependent survival by lowering competition for algae.Together, these results suggest that a structural change in the validated daphnia model, i.e., including predation, and not just a change in parameter values is required to reproduce population oscillations observed in laboratory and natural settings.This finding is important in the context of our previous and current ongoing efforts, since oscillations were not observed in our previous daphnia population experiments [5].Although the experiments performed in this work did not involve toxins, the species investigated are commonly used in toxicity assessments.Thus, our results can be used to provide a baseline to compare effects in a toxicity setting in future work.

Figure 1 .
Figure 1.Plots of forward solutions for the Logistic curve for the three replicates of the data.Replicate one is on top and three is on the bottom.The lighter and darker shades of grey represent the 95% and 68% confidence bars on the model solution, respectively.The algae population is represented as cells × 10 7 /L.Data points are shown as "*".

Figure 2 .
Figure 2. Plots of forward solutions for the discrete Euler-method logistic (DEL) please confirm.model for the three replicates of the data from left to right.Replicate one is on top and three is on the bottom.The lighter and darker shades of grey represent the 95% and 68% confidence bars on the model solution, respectively.The algae population is represented as cells × 10 7 /L.Data points are shown as "*".

Figure 3 .
Figure 3. Plots of forward solutions for the Gompertz curve for the three replicates of the data.Replicate one is on top and three is on the bottom.The lighter and darker shades of grey represent the 95% and 68% confidence bars on the model solution, respectively.The algae population is represented as cells × 10 7 /L.Data points are shown as "*".

Figure 4 .
Figure 4. Plots of forward solutions for the Bernoulli curve for the three replicates of the data.Replicate one is on top and three is on the bottom.The lighter and darker shades of grey represent the 95% and 68% confidence bars on the model solution, respectively.The algae population is represented as cells × 10 7 /L.Data points are shown as "*".

Figure 5 .
Figure 5. Simulations of the coupled daphnia and green algae dynamics model resulting in steady state behavior (Left, c = 0.01) and sustained oscillations (Right, c = 0.04).
) where, for our own examples, n = 23 and p = 3 or 4 is the number of model parameters.Both θn

Table 1 .
Corrected Akaike Information Criterion scores for each model and replicate.

Table 2 .
Model comparison ûn scores for the continuous Logistic model compared to Bernoulli model for each replicate.We also chose to fix the initial condition X 0 at the seed population value to enunciate model improvement if X 0 was unalterable.Note that values less than 2.706 indicate the restricted model is better.

Table 3 .
β estimate and standard error for the Bernoulli model.
1/P t if i ≤ 4. In this model, we changed the functional form of the Daphnia model based on the assumption that the negative density-dependent fecundity and survivorship effects that daphnids experience are driven by a lack of food in the form of algae, represented as 1 P t and its time-delayed analogue 1 P t−τ .The Daphnia matrix model is otherwise unchanged.

Table 4 .
Summary of Daphnia magna and algae model parameters and variables.