The effects of regular seasonal variation on population models of two or three interacting species have been studied extensively (Kot & Schaffer 1984; Doveri et al. 1993; Rinaldi et al. 1993; Steffen et al. 1997; Huppert et al. 2005). These studies have shown that periodically forced populations can display a rich repertoire of dynamical behaviours, including simple and complex periodic cycles, quasi-periodicity and chaos (Rinaldi et al. 1993; King & Schaffer 1999; Vandermeer et al. 2001). However, the parameter space in which chaotic behaviour occurs is usually small. Typically, the population dynamics show repeatable patterns. Slow-growing species may fluctuate on seasonal or multi-annual time scales, as exemplified by the famous cycles of voles and lemmings at northern latitudes (Stenseth 1999; Turchin 2003). Fast-growing species such as bacteria or plankton may display multiple ups and downs per year. The frequency of the population fluctuations can be remarkably persistent as a result of frequency locking (e.g. Scheffer et al. 1997; Vandermeer et al. 2001). Seasonal forcing tends to ‘lock’ the frequency of population oscillations, such that populations oscillate at the same frequency as the seasonal cycle or integer multiples of it.

While many theoretical studies have examined effects of seasonality on model systems of only a few species, seasonal forcing of multi-species communities has received surprisingly little theoretical attention (Ebenhöh 1992). Yet bacterial, plankton and insect communities may contain tens, hundreds and sometimes even thousands of species (Hutchinson 1961; Erwin 1982; Irigoien et al. 2004; Venter et al. 2004). Generally speaking, multi-species models display more complex dynamics than models with only two or three species (May 1973; Ellner & Turchin 1995; Huisman & Weissing 2001). From a conceptual perspective, multi-species food webs can be interpreted as systems with several interacting oscillations (e.g. several predator–prey cycles). Coupled oscillations are known to generate complex dynamics, including chaos (Hastings & Powell 1991; Vandermeer 1993, 2004; Huisman & Weissing 2001; Benincà et al. 2008). The prevalence of complex dynamics is of interest, because these non-equilibrium dynamics may help to sustain the biodiversity of natural communities (Armstrong & McGehee 1980; Huisman & Weissing 1999; Brose 2008) and also because complex dynamics can induce regime shifts in ecosystems, with important implications for their management (Scheffer et al. 2001; Hsieh et al. 2005; Ives et al. 2008).

These dynamics can be illustrated by Poincaré maps sampling the model communities once per year for many consecutive years. Model communities with a periodicity of 1 year return to exactly the same species composition year after year, which appears as a single point on the Poincaré map. Communities with a periodicity of N years produce N points on the Poincaré map, quasi-periodicity produces a closed curve (figure 3a), while chaos produces a complex fractal structure (figure 3b).

Many of the model communities exposed to seasonal forcing displayed chaos with remarkable synchronization patterns at the species level (figure 4). The species fluctuations are irregular, yet these irregular fluctuations are squeezed within the seasonal cycle. As a consequence, species enter the winter season in different proportions, and this affects the species composition of the next spring bloom. For instance, figure 4a shows a typical phytoplankton spring species. It reaches peak abundance in March, although its peak abundance varies from year to year, and some years it does not peak in spring at all. Figure 4b shows another phytoplankton species from the same plankton community. This species could be called a typical summer species. It is present every summer. However, some years it peaks twice, with a first peak in May–June and a second smaller peak in September. In other years, it peaks in September only. The zooplankton species show similar seasonal patterning. For instance, some zooplankton species are mainly present in winter (figure 4c), while others dominate during the summer period (figure 4d). The example in figure 4d is particularly interesting. In some years, this zooplankton species shows little variability from March to September, while in other years, it fluctuates wildly during the same period. Accordingly, the species composition in our model communities shows distinct patterns of seasonal organization, but with strong year-to-year variability.

Which species traits and environmental conditions are responsible for the widespread chaotic dynamics in our model communities? A complete answer to this question is beyond the scope of this paper. However, some insight can be obtained by modifying the model assumptions systematically. This shows that more than 50 per cent of the model simulations produced chaos when using our default parameter settings (table 2, first row). The occurrence of chaos was not very sensitive to the relative magnitude of intraspecific versus interspecific phytoplankton competition (table 2). In contrast, modifying zooplankton predation had a striking effect on the occurrence of chaos. When zooplankton was removed from the model, very few simulations showed chaotic dynamics and they did so only under seasonal forcing (table 2). Similarly, inefficient zooplankton grazing and specialist zooplankton reduced the occurrence of chaos. This shows that predator–prey oscillations, and the nature of predation, played a key role in the generation of complex dynamics in our model communities.

Productivity also had a clear effect on the occurrence of chaotic dynamics. At low productivity (K = 2 mg l⁻¹), stationary dynamics prevailed in constant environments, simple periodic dynamics prevailed in seasonal environments and chaos occurred only in a few model communities with strong seasonal forcing (table 2; see also appendix S1 in the electronic supplementary material). Chaos was widespread at intermediate productivity (K = 5 and K = 10 mg l⁻¹). At high productivity (K = 20 and K = 50 mg l⁻¹), the occurrence of chaos declined and the population dynamics often shifted to periodic cycles in both constant and seasonal environments.

In all cases summarized in table 2, seasonal forcing increased the occurrence of chaos. To investigate this aspect in further detail, we estimated whether the amplitude of seasonal forcing affected the occurrence of chaos in our model communities (figure 5). We focused on the intermediate productivities (K = 5 and K = 10 mg l⁻¹). At K = 5 mg l⁻¹, the amplitude of seasonal forcing increased the occurrence of chaos (figure 5a; linear regression: R² = 0.45, N = 11, p = 0.024). At K = 10 mg l⁻¹, mild forcing (0.1 < a < 0.4) caused a slight increase in the probability of chaos, but when the amplitude of seasonal forcing was further increased (a > 0.6), the probability of chaos declined (figure 5b; quadratic regression: R² = 0.66, N = 11, p = 0.013). We further explored the predictability of these communities by calculating their Lyapunov exponents. A positive Lyapunov exponent indicates chaos. The inverse value of the Lyapunov exponent is often used as a simple metric of the predictability of chaotic systems (Strogatz 1994). In those simulations that displayed chaotic dynamics, the magnitude of the Lyapunov exponent was not affected by the amplitude of seasonal forcing (figure 5c,d). This indicates that the predictability of the chaotic plankton communities was neither enhanced nor reduced by a stronger seasonality. However, the median values of the Lyapunov exponents were significantly higher at K = 10 mg l⁻¹ than at K = 5 mg l⁻¹ (figure 5c,d; t-test: t = −3.77, d.f. = 20, p < 0.002), which indicates that the predictability of the model communities was affected by productivity.