Stochastic dynamics of Cyanobacteria in long‐term high‐frequency observations of a eutrophic lake

Concentrations of phycocyanin, a pigment of Cyanobacteria, were measured at 1-min intervals during the icefree seasons of 2008–2018 by automated sensors suspended from a buoy at a central station in Lake Mendota, Wisconsin, U.S.A. In each year, stochastic-dynamic models fitted to time series of log-transformed phycocyanin concentration revealed two alternative stable states and random factors that were much larger than the difference between the alternate stable states. Transitions between low and high states were abrupt and apparently driven by stochasticity. Variation in annual magnitudes of the alternate states and the stochastic factors were not correlated with annual phosphorus input to the lake. At daily time scales, however, phycocyanin concentration was correlated with phosphorus input, precipitation, and wind velocity for time lags of 1–15 d. Multiple years of high-frequency data were needed to discern these patterns in the noise-dominated dynamics of Cyanobacteria. Blooms of harmful Cyanobacteria have adverse effects on water quality, survival of aquatic animals, and human health in lakes, reservoirs, and coastal oceans (Carmichael and Boyer 2016; Huisman et al. 2018). Blooms are associated with high inputs of nutrients interacting with climate change (Michalak et al. 2013). Growing demand for food and a shift toward more *Correspondence: steve.carpenter@wisc.edu Associate editor: Monika Winder Author Contribution Statement: The study was conceived by S.R.C., E.H.S., and P.C.H., the automated sensor system was designed by P.C.H., modeling was conceived and conducted by S.R.C., B.M.S.A., M.S., and E.V.N., and the paper was written by all authors. Data Availability Statement: Data are available in the North Temperate Lakes Long-Term Ecological Research database (http://lter.limnology. wisc.edu). Additional Supporting Information may be found in the online version of this article. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

meat consumption by a growing human population are increasing the exports of phosphorus from agriculture to surface waters (MacDonald et al. 2011). At the same time, precipitation is becoming concentrated into larger storms, thereby increasing erosion and nutrient transport to surface waters (Carpenter et al. 2018). Cyanobacterial dominance is further promoted by rising temperatures (Johnk et al. 2008;Kosten et al. 2012).
Within a eutrophic lake or reservoir, day-to-day fluctuations of Cyanobacteria can be highly variable (Scheffer et al. 2003;Huisman et al. 2018). It is these short-term fluctuations that affect daily human use of the water supply for drinking, recreation, or harvesting fish. Models of phytoplankton blooms in the context of nutrient and food web dynamics show diverse behaviors including multiple stable points and cycles, depending on the current rate of nutrient input and state of the consumer food web (Scheffer et al. 2000). This dynamical diversity combined with fluctuations of the physical environment suggests that planktonic ecosystems may never be at equilibrium, although interannual and multilake landscape patterns are predictable (Padisák et al. 1993;Scheffer et al. 2003).
To compare the roles of equilibria vs. environmental variability for Cyanobacteria dynamics in surface waters, we studied fluctuations of phycocyanin, a characteristic pigment of Cyanobacteria, at 1-min intervals using sensors suspended from a centrally located buoy during the ice-free seasons of 11 yr in eutrophic Lake Mendota, Wisconsin. The time series were analyzed using drift-diffusion-jump models (Carpenter and Brock 2011) (Supporting Information). In each year, we found abrupt shifts between stable regimes of low and high phycocyanin. However, stochasticity at the scale of minutes was larger than the gap between alternate states, and sometimes linked to rapid shifts between regimes.

Lake Mendota
Limnological studies of Lake Mendota, Wisconsin, U.S.A. (43 0.6 0 N, 89 25 0 W, 3.98 × 10 7 m 2 surface area, 12.3 m mean depth) began in 1875 and the lake contained bloom-forming Cyanobacteria by the 1880s (Carpenter et al. 2006). Eutrophication is driven primarily by runoff of phosphorus from agriculture, with some contributions from urban runoff (Lathrop et al. 1996;Lathrop and Carpenter 2014). Composition and seasonality of the bloom-forming Cyanobacteria taxa were described by Lathrop and Carpenter (1992).

Automated sensor
Concentrations of phycocyanin and chlorophyll a were monitored every minute at a central station in Lake Mendota during the ice-free seasons of 2008-2018. Data are available in the North Temperate Lakes Long-Term Ecological Research program database (https://lter.limnology.wisc.edu). Phycocyanin and chlorophyll data were recorded every minute using Turner Designs Cyclops 7 sensors suspended 1.0 m below the buoy.
The duration of the measurements is around 175 d each year, offering rather long time series of approximately 250,000 observations annually. Phycocyanin sensor readings were expressed as relative fluorescence units (RFU), which are directly related to concentrations of the respective pigments (Pace et al. 2017).

Time-series analysis
We fit a drift-diffusion-jump model (Johannes 2004;Carpenter and Brock 2011;Anvari et al. 2016) to summarize the dynamic patterns of phycocyanin observations over 1-min intervals of each ice-free season. For each ice-free season, the sequence of analysis was (1) convert log 10 (phycocyanin) to z-scores (y z,t = y t − y À Á =s where y is the sample mean and s is the sample standard deviation), (2) from the z-scores calculate standardized levels of phycocyanin, (3) fit a drift-diffusion jump model to the standardized levels, and (4)   selected by Akaike Information Criterion (AIC), was 1 in all cases. As a smoothed indicator of phycocyanin, we used the time-varying intercept divided by its standard error b t /s t . We refer to the time series b t /s t as the "standardized level of phycocyanin." We used the standardized levels of phycocyanin as input (x t ) to fit the drift-diffusion-jump model (Carpenter and Brock 2011) (Supporting Information). D 1 (x) is the deterministic core of the dynamics as a function of x called the "drift" in stochastic dynamic modeling. Its roots D 1 (x) = 0 are the equilibria. D 2 (x) is the Gaussian component of noise expressed as a standard deviation that is function of x, known as the "diffusion." J(x) is the Poisson component of noise known as the "jump." J(x) is characterized by the Poisson jump rate λ(x) and jump size ξ which is assumed to be normally distributed as ξ~N(0, σ ξ 2 ) where σ ξ 2 is called jump amplitude. We estimated drift, diffusion, and jump functions, that is, D 1 (x), D 2 (x), jump rate λ(x), and jump amplitude σ ξ 2 , by nonparametric regression using the algorithm of Johannes (2004) implemented in R (Carpenter and Brock 2011). Equilibria and the contribution of jumps to total variance (diffusion + jumps) were computed and  back-transformed to the original units, log 10 (phycocyanin), as shown in Supporting Information.

Results
Observed time series were highly variable (Fig. 1A), spanning a range of about 20 times the standard deviation within each year. The standardized level (time-varying intercept divided by its standard error) had clearly-defined intervals of low or high values punctuated by sharp shifts (Fig. 1B). High values of b t /s t result from high fitted intercepts b t with small standard errors s t during blooms.
The drift function for 2017 crossed the zero line three times, indicating three equilibria: a stable equilibrium at low levels of phycocyanin, an unstable equilibrium at intermediate levels, and another stable equilibrium at high levels of phycocyanin ( Fig. 2A). The diffusion function, or Gaussian component of Eq. 1, was relatively large, consistent with the high variability apparent in the time series (Fig. 2B). Figure 2C presents the equilibrium values on a plot of the standardized level vs. day of the year for 2017. There are frequent shifts between low and high phycocyanin states, showing the abrupt onset and collapse of blooms.
Values of the equilibria (log 10 of phycocyanin RFU) varied about 10-fold from year to year (Fig. 3A). Diffusion (in the same units) was substantially higher around the unstable equilibria than around the stable equilibria (Fig. 3B). Moreover, diffusion varied a million-fold among years, especially around the unstable equilibria. The percentage of total variance due to the Poisson component, or jump function (Supporting Information), also varied among years, ranging from near 0% to about 80% (Fig. 3C). The role of jumps was consistently greater around the unstable equilibria.
No significant correlations were found between annual phosphorus load (Supporting Information Fig. S7) and any of the annual parameters of phycocyanin plotted in Fig. 3. However, the sample size was small (11 yr) and these correlations have low power. At daily resolution, sample size is higher and there are significant correlations of phycocyanin time series with phosphorus load >9 d earlier, precipitation 12-16 d earlier, and wind velocity 1-5 d earlier (Supporting Information Fig. S1).
The range of outcomes is illustrated by the contrast between 2010 (Fig. 4A) and 2016 (Fig. 4B). In 2010, the equilibria were relatively low, diffusion was relatively high, and jumps were dominant (Fig. 3). In 2016, the equilibria were relatively high and diffusion was relatively low (Fig. 3).

Discussion
Phycocyanin dynamics in Lake Mendota were dominated by stochasticity. Two stable equilibria were found in each year, but the gap between the stable equilibria was small relative to the noise. The differences among years are larger than the differences between stable equilibria seen in any individual year. The magnitude of stochasticity also showed large differences among years.
The interannual variation of equilibria was not correlated significantly with phosphorus loading, a major driver of Cyanobacteria in Lake Mendota (Lathrop and Carpenter 2014). Analyses of a longer time series (21 yr) of phosphorus loading and manual counts of Cyanobacteria quantified the relationship of phosphorus inputs to Cyanobacteria mean biomass and bloom frequency (Stow et al. 1997;Lathrop et al. 1998). Daily sampling during the ice-free season of 1993 showed that blooms of Cyanobacteria followed changes in meteorological conditions and grazer biomass with time lags of 0-9 d (Soranno 1997). In this study, we found similar effects of phosphorus load, wind velocity, and precipitation on phycocyanin concentrations at lags of 5-14 d (Supporting Information Fig. S1). Thus some of the variability of Cyanobacteria is related to storm-driven pulses of phosphorus input, wind-driven mixing, and grazing.
The paradox of enrichment suggests that Cyanobacteria dynamics exhibit cycles of increasing magnitude as nutrient supply increases (Rosenzweig 1971). More comprehensive models show the possibility of multiple stable states or cycles, depending on the interactions of nutrient input, grazers, and fish (Scheffer et al. 2000). Syntheses of Cyanobacteria dynamics in many lakes attribute their variability to weather-induced fluctuations, spatial complexity of the turbulent mixed-layer environment, and endogenous cycles or chaos (Padisák et al. 1993;Scheffer et al. 2003).Our data indicate two stable states of phycocyanin concentration that vary among years. During summer stratification, average Secchi disk transparency ranges from about 1.5 to 3 m (Lathrop et al. 1996), depending on phosphorus inputs and density of large-bodied grazers especially Daphnia pulicaria (Lathrop et al. 1999). At such low transparency, much of the mixed layer is dark and phytoplankton growth may be limited by low irradiance. Therefore the high stable state of Cyanobacteria concentration may be determined by a balance of growth and self-shading. The causes of the low stable state may be related to a balance of growth and grazing losses.
!The data we report were collected at a single station in the lake and therefore cannot assess the spatial variability of Cyanobacteria in surface waters. Studies with boat-mounted sensors in Lake Mendota show that spatial variability is large (Loken et al. 2019). We cannot determine the contribution of spatial variation to the stochastic models that we fit to these time series. Better integration of spatial dynamics with long-term variability is an important need for future studies of Cyanobacteria.
This decade of high frequency time series observations reveals a pattern of low and high stable states of Cyanobacteria concentration within a regime of high stochasticity. A short-term study or low-frequency sampling would miss this overarching pattern of Cyanobacteria dynamics. By combining high frequency and long-term observations with driftdiffusion-jump models, this study has revealed fundamental regularities in the stochastic dynamics of Cyanobacteria.