Journal of Theoretical Biology

In food web models that include more than one prey type for a single predator, it is common for the predator’s functional response to include some form of switching—preferential consumption of more abundant prey types. Predator switching promotes coexistence among competing prey types and increases diversity in the prey community. Here, we show how the dynamics of a diamond-shaped food web model of a marine plankton community are sensitive to a parameter that sets the strength of predator switching. Stronger switching destabilizes the model’s coexistence equilibrium and leads to the appearance of limit cycles. Stronger switching also increases the evenness of the asymptotic prey community and promotes synchrony in the dynamics of disparate prey types. Given the dependence of model behavior on the strength of predator switching, it is important that modelers carefully consider the parameterization of functional responses that include switching.


Introduction
In mathematical models of ecological community dynamics, a predator's functional response (Solomon, 1949;Holling, 1959a,b) describes the predator's per capita grazing rate as a function of prey density. The choice of functional response plays a crucial role in determining the asymptotic dynamics of a model, including the equilibrium abundance and stability of predator and prey populations. Reflecting their importance, functional responses for many kinds of consumers have been the object of intense study in ecology. Uiterwaal et al. (2022) documented over 2000 instances of functional responses that have been estimated in the scientific literature. Methods for estimating functional responses from data continue to be developed (Coblentz and DeLong, 2021), and mathematical analyses reveal that subtle differences in the form of the functional response can have out-sized effects on community dynamics (Fussmann and Blasius, 2005).
When predators consume more than one type of prey, the choice of predator functional response is even more complex because the function must account for the predator's choice between multiple potential prey. This complexity is especially relevant, for example, in models of marine planktonic food webs, which often include many coexisting phytoplankton types. The number of different functional responses available to modelers is large, and choosing between them can be difficult since the assumptions implicit in each may not be clearly stated and the dynamics of the models themselves are often sensitive to the details of the functional response (Fussmann and Blasius, 2005). Ideally, * Corresponding author at: University of California Santa Barbara, Santa Barbara, CA, USA.
E-mail address: karchibald@ucsb.edu (K.M. Archibald). the choice of a functional response should follow careful consideration of the ecological phenomena that are relevant to the question at hand, as well as the mathematical characteristics of the functional response and how they may impact the model's dynamics.
One phenomenon that is often included in models of marine food webs with multiple phytoplankton types is switching, which occurs when a zooplankter's grazing rate on a phytoplankton type is disproportionate to that type's frequency in the environment (Gentleman et al., 2003). Switching represents a zooplankter's ability to preferentially consume more abundant phytoplankton types. Switching is frequently observed in laboratory and field studies in marine ecosystems (Murdoch, 1969;Murdoch et al., 1975;Hughes and Croy, 1993;Kiørboe et al., 1996;Kempf et al., 2008) and modeling work has shown that switching can stabilize dynamics and increase diversity (Armstrong, 1999;Morozov, 2010;Prowe et al., 2012;Vallina et al., 2014). In this way, predator behavior can influence the asymptotic dynamics of plankton communities and interact with other factors like competition for shared resources (reviewed by Tilman (1982)), fluctuations in abiotic conditions (Levins, 1979;Ebenhöh, 1988), and the mathematical form of transfer functions and closure terms (Sjöberg, 1977;Steele and Henderson, 1992;Ruan, 2001).
A number of early multi-prey functional responses that included switching suffered from an undesirable mathematical characteristic, called antagonistic grazing (Tilman, 1982;Vallina et al., 2014). Antagonistic grazing occurs when a zooplankter's total consumption rate is https://doi.org/10.1016/j.jtbi.2023.111536 Received 12 January 2023; Received in revised form 27 April 2023; Accepted 11 May 2023 higher when feeding on a single phytoplankton prey than it would be if the zooplankter were feeding on the same total biomass divided among several phytoplankton types. To fix the antagonistic grazing problem, Vallina et al. (2014) developed the so-called ''Kill-the-Winner'' (KTW) functional response to describe switching in a distributed grazing model. In the KTW response, the total grazing rate is a function of the total prey biomass, ensuring that the total grazing rate is independent of the number of phytoplankton types or their distribution in the environment. The total grazing pressure is then distributed across different phytoplankton types based on the relative abundance of each type. Since its publication, KTW has become a popular way of including switching in many models of a variety of systems (Guyennon et al., 2015;Egilmez and Morozov, 2016;Ward and Follows, 2016;Baudrot et al., 2016;Smith et al., 2016;Cadier et al., 2017;Vallina et al., 2017;Tanioka and Matsumoto, 2018;Nissen et al., 2018;Chen et al., 2019).
In a community with phytoplankton types with densities , the KTW functional response takes the form where is the per capita consumption rate of the zooplankter on phytoplankton type (Vallina et al., 2014). The model posits that grazing is a saturating function of total prey abundance (with a halfsaturation prey abundance, ), such that zooplankton follow a Holling Type II functional response (Holling, 1959a). The KTW functional response also accommodates both prey preference and prey switching behaviors: , is the grazer's preference (i.e., not frequency dependent) for phytoplankton type , and governs the steepness in transition of the zooplankter's grazing on the most abundant phytoplankton type (i.e., ''switching''; Fig. 1). When = 1, no switching occurs, but when ≫ 1, the zooplankter exhibits strong preferential grazing on the most abundant phytoplankton type. 1 Finally, the parameter max is the zooplankter's maximum grazing rate, which is achieved when at least one of the phytoplankton types is extremely abundant. The totality of the zooplankter's grazing is distributed across all phytoplankton types according to their frequency, the zooplankter's preferences, and the strength of switching. Because our focus here is on the effects of predator switching (governed by the value of ), we set = 1 for all in the analyses that follow to eliminate the effects of zooplankton preference (but see Supplemental Fig. S1 for the effects of preferences on the shape of the KTW functional response).
The key characteristic that distinguishes KTW from other multi-prey functional responses is how the zooplankter's total ingestion rate is calculated (Vallina et al., 2014). In KTW, total ingestion rate depends only on total prey abundance and not on the distribution of phytoplankton types. This ensures the functional response does not exhibit antagonistic grazing. Switching, therefore, is implemented as a separate term that distributes the total grazing pressure across all phytoplankton types based on their frequency. This frequency-dependent grazing represents a behavioral change in the grazer (e.g. feeding strategy, capture efficiency) in response to variance in the prey community (Murdoch, 1973;Kiørboe et al., 1996).
The KTW functional response is phenomenological; it is intended to account for a particular phenomenon, rather than being derived from specific biological mechanisms. This includes the parameter , which describes the strength of predator switching in the model (Vallina et al., 2014). Due to the difficulty of directly measuring , there is a lack of data in the oceanographic literature describing the strength of switching or the variability of between different biogeochemical regimes or for different community assemblages. Most studies that utilize KTW have defaulted to a value of = 2 based on precedent (Vallina et al., 2014). In cases where the effects of switching are evaluated, it is most common to see switching treated as a binary characteristic; that is, comparing model behavior when there is no switching to a case when there is some (arbitrary) degree of switching. One notable exception is Smith et al. (2016), who have shown that increasing results in higher phytoplankton diversity in size-structured nutrient-phytoplankton-zooplankton-detritus (NPZD) models.
Here, we investigate how the strength of predator switching affects the asymptotic dynamics of a plankton food web model with the KTW functional response. We focus on a broad range of parameter values to explore how both the stability of the model equilibria and the characteristics of the phytoplankton community (e.g., evenness, synchrony) change as the strength of predator switching increases. We find that stronger switching destabilizes the model and leads to the appearance of limit cycles. Stronger switching also results in increased evenness in the phytoplankton community and promotes synchrony in the temporal variability of different phytoplankton types.

The model
We study a three-level food web in which phytoplankton types with population sizes ( ∈ 1, … , ) compete for a shared essential nutrient resource and are exposed to predation from a single zooplankton predator . When = 2, this model produces a ''diamond'' food web (Evans, 1988).
We model these competition and predation dynamics as taking place within a chemostat with a dilution rate . The nutrient flows into the chemostat at a concentration 0 , where it may be taken up by phytoplankton or diluted out of the system. We assume that phytoplankton take up the nutrient resource following a saturating functional response (Monod, 1949) with a specific uptake rate and half-saturation constant . The dynamics of the nutrient resource are given by: (2) K.M. Archibald et al.

Fig. 2.
The per-capita grazing rate on phytoplankton type 1 ( 1 (1, 2 )) as a function of the density of phytoplankton type 2 ( 2 ), when 1 is held constant ( 1 = 1) for selected values of . The functional response exhibits synergistic grazing because 1 ( 1 , 2 )∕ 2 > 0 for some 2 >= 0; this only occurs when > 1. Phytoplankton dynamics are governed by the balance between nutrient uptake, loss through dilution out of the chemostat, and mortality due to predation at a species-specific grazing rate : with given by Eq.
(1). Note that we are measuring the biomass of all living organisms ( and ) in units of nutrients. Conversion efficiencies both here and in zooplankton grazing terms are neglected.
Finally, zooplankton growth arises from the balance of grazing on phytoplankton and losses due to dilution:

Dimensionless equations
In order to simplify our analysis, we rescaled the model's variables as follows: Substitution yields the dimensionless model: with the rescaled KTW functional response and rescaled parameters Because our model components are contained in a chemostat, conservation of mass dictates that, asymptotically, K.M. Archibald et al. With (8), we can reduce our system of equations to + 1 equations: In cases where > 1, we established a hierarchy among the phytoplankton types so that type is competitively superior to type when < . We achieve this hierarchy by linearly scaling the nutrient uptake rate as in Table 1 (so that 1 > 2 > ⋯ > ). All phytoplankton types were assigned the same nutrient half saturation constant ( = 0.01; Table 1). We can therefore measure competitive strength as the difference between the values of two phytoplankton. The phytoplankter with the highest value ( 1 ) is competitively dominant because, in the absence of predation, it can draw the nutrient level in the chemostat down to the lowest level at equilibrium (Tilman, 1977) (Fig. S2).
We analyze model (9) by first considering the two-phytoplankton case ( = 2) in order to understand how the strength of zooplankton switching ( ) impacts the dynamics of phytoplankton coexistence at steady-state. We then ask how the model's asymptotic behavior depends upon , phytoplankton traits (especially , the prey resource uptake rate), and phytoplankton diversity. Throughout the analysis, we make note of two characteristics of phytoplankton community dynamicsevenness and synchrony. Community evenness ( ′ ) is defined as the ratio of the simulated diversity to the maximum possible diversity in a community with the same number of phytoplankton types (Pielou, 1966), where is the biomass ratio of phytoplankton type to the total phytoplankton community (i.e., = ∕ ∑ ). Community synchrony, on the other hand, describes the degree to which the temporal variance of different plankton types are correlated. For example, in a highly synchronous community, all phytoplankton types will reach maximum or minimum biomass at the same time. In an asynchronous community, different types will reach their peak biomass at different times, possibly maintaining a constant community-level total biomass as a result of compensatory biomass variance. Synchrony between phytoplankton types determines whether the variance of ecosystemlevel ecological functions (e.g. productivity) reflects the variance at the level of individual phytoplankton types.

KTW promotes coexistence between prey species
The presence of switching in the KTW functional response allows two competing phytoplankton types to coexist in situations that would otherwise lead to competitive exclusion. Consider that phytoplankton compete directly with one another for resources (in our model, a shared nutrient pool). Therefore, in the absence of predation, the weaker competitor would be excluded from the system following classic R * theory (Tilman, 1977). Because they are grazed by a shared zooplankton predator, phytoplankton are also subjected to indirect competition (Holt, 1977), wherein the presence of a second phytoplankton prey supports a numerical increase in the zooplankton, which in turn increases predation pressure on the first phytoplankton type. By focusing grazing pressure on numerically dominant phytoplankton types, switching relieves the effect of indirect competition on weaker phytoplankton competitors while simultaneously intensifying the effect on stronger competitors. We refer to this phenomenon as synergistic grazing -a positive relationship between the abundance of a focal prey type and the per-predator consumption rate of an alternative prey type when the abundance of the alternate type is low (Fig. 2). In other words, a numerical increase in one phytoplankton type leads to greater per-zooplankter grazing on an alternative phytoplankton type. A functional response is said to exhibit synergistic grazing if ( , )∕ for some >= 0. This somewhat counter-intuitive model behavior arises from the fact that the KTW function is designed to maintain constant grazing rates as a function of total prey biomass. When a second phytoplankton type is introduced to the system at a lower abundance than the first, the KTW function generates additional predation on the first phytoplankton type, especially when is large. This creates a refuge from predation at low abundances, allowing multiple phytoplankton types to coexist, so long as switching is present ( Fig. 3b-d).
The role of switching in promoting coexistence is clearly demonstrated by examining the hypothetical scenario in which a phytoplankton type is trying to invade a system in which its competitor is already present. Consider the equilibrium where , 1 , and are at their steadystate values ( * , * 1 , * ). The invasion growth rate ( 0 ) of 2 at this equilibrium point can be derived from Eq. (9a) (see Appendix A): If 0 > 0, then 2 is able to invade the equilibrium and coexistence is possible. Since 0 is strictly larger when switching is included in the functional response ( > 1), coexistence will be possible for a larger range of parameter values when the zooplankter grazer exhibits prey switching. Interestingly, the simple presence or absence of switching influences the coexistence criteria, but the magnitude of does not, as the equilibrium values * , * , and * do not depend on in Eq. (11). So while, as we will demonstrate next, the strength of switching has a strong effect on the evenness of the steady-state prey community, it has no effect on the number of coexisting phytoplankton types. Indeed, as long as some degree of switching is present in the functional response, an arbitrary number of phytoplankton types can coexist at equilibrium at sufficiently low abundances (provided that those types are capable of positive net growth in isolation).

Strong switching is destabilizing.
Model (9) exhibits two kinds of asymptotic behavior: stable equilibria and limit cycles (Fig. 4). We did not find bistability in our model; for each combination of parameters, only one type of asymptotic behavior was observed. Stronger zooplankton switching destabilizes stable equilibria in favor of limit cycles (Fig. 5a). Limit cycles also emerge when the number of coexisting phytoplankton species is high (Fig. 6). This destabilization appears to arise from the intensification and synchronization of predation pressure. For larger values of , the predator's switching function becomes increasingly steep, leading to sharper transitions in predation intensity between numerically dominant and sub-dominant phytoplankton types. This mechanism functions to homogenize phytoplankton population sizes by disproportionately punishing numerically dominant phytoplankton types. Similarly, when the phytoplankton community is subdivided into many different types, their population sizes are more similar. This similarity causes the phytoplankton community to respond synchronously to changes in the zooplankton population, which can act to amplify the system's propensity for predator-prey cycles. This hypothesis is further supported by the appearance of limit cycles when the phytoplankton's growth rates are more similar (i.e., phytoplankton populations with smaller competitive differences tend to exhibit limit cycles and populations with larger competitive differences tend to exhibit stable equilibria; Fig. 5b).

Switching promotes evenness and synchrony in prey populations.
Stronger switching also promotes evenness in the prey community (Fig. 7). The higher the degree of switching, the more the zooplankton preferentially consume the numerically dominant phytoplankton types. This top-down control acts to suppress the most numerically dominant phytoplankton, alleviating competition with less dominant phytoplankton types, whose numbers subsequently increase. As a consequence, the sub-dominant competitors have population peaks that lag slightly Fig. 5. Asymptotic behavior of 1 and 2 as a function of (a) and (b) the competitive difference between 1 and 2 ( 1 − ). The appearance of limit cycles is associated with stronger switching and reduced competitive differences between phytoplankton types. In (a) 1 = 12 and 2 = 4; in (b) 2 = 4 and = 2. behind the peak of the dominant competitor ( Fig. 3c-d). The ability for switching to promote evenness is quite high; community evenness increases rapidly for values of between 1 and 2, with evenness converging to the theoretical maximum value of 1 for very strong switching (Fig. 7). This increase in evenness occurs independently of the total prey community biomass, which does not depend on (Fig. 7).
Interestingly, the presence of switching in the functional response appears to encourage synchronous dynamics in phytoplankton prey populations (Fig. 3d). Populations that are out of phase maximize the numerical differences in concentration since the maximum of one type coincides with the minimum of another type. Switching, through the preferential consumption of the numerically dominant phytoplankton type at any given point in the limit cycle, acts to counter these imbalances and synchronize prey population dynamics. The variance of the total phytoplankton community (summed across all types) is high if population fluctuations are synchronous (i.e. in phase), as opposed to asynchronous (out of phase) fluctuations that may compensate for variability at the population level and maintain a fixed community-level biomass. As a result, switching does not function to stabilize overall ecosystem function (e.g., productivity, biomass) at the community level.

Discussion
Grazing functional responses that include switching, such as KTW, are used in many types of ecological models to promote diversity in the prey community (Guyennon et al., 2015;Egilmez and Morozov, 2016;Ward and Follows, 2016;Baudrot et al., 2016;Smith et al., 2016;Cadier et al., 2017;Vallina et al., 2017;Tanioka and Matsumoto, 2018;Nissen et al., 2018;Chen et al., 2019). Our analysis of a marine food web model has revealed that in addition to increasing diversity (and evenness) in the phytoplankton community, switching strength can also change the asymptotic behavior of food web models. Stronger switching promotes instability of the model's equilibrium, leading to the emergence of limit cycles. These limit cycles are particularly important since switching appears to drive synchrony between prey types, thereby reducing the stability of ecosystem provisioning (e.g., productivity or biomass).
Interestingly, we also observed that higher prey diversity was associated with the appearance of limit cycles in the model. This negative relationship is perhaps counter-intuitive in the context of traditional views of diversity and stability, in which increased diversity promotes increased stability of ecosystem function (Tilman et al., 1998;Elton, K.M. Archibald et al. Fig. 7. Expected evenness (gray) and total phytoplankton biomass (green) in the asymptotic prey community as a function of switching strength ( ). Shaded regions indicate the range of values in cases where the asymptotic behavior of the model is limit cycles. Here, 1 = 16, 2 = 4, = 2, and = 16.
1958; MacArthur, 1955). Diversity is often assumed to provide ecological redundancies that decrease the volatility of the ecosystem. In this case, however, the synchrony between the phytoplankton types means that there is no compensation between different phytoplankton types. And since increasing the number of prey types can result in the appearance of limit cycles, increased diversity drives increased variability at both the population and community levels.
These results provide important context for the ecological role of predator switching, which is frequently observed in both laboratory and field studies (Murdoch, 1969;Murdoch et al., 1975;Hughes and Croy, 1993;Kiørboe et al., 1996;Kempf et al., 2008). As an ecological phenomenon, switching may arise via multiple mechanisms that all have the same emergent behavior: frequency-dependent differential grazing on multiple available prey types. The classical view of switching is that it arises from predator behavior, such as changes in feeding strategy or time budget optimization in response to a variable prey community (Paffenhöfer, 1984;Strom et al., 2000;Mariani and Visser, 2010). Alternatively, the zooplankton represented in this model could be viewed as a generalized micro-predator community rather than a single species. Under this scenario, the switching ''behaviors'' represent internal changes in the composition of the zooplankton community as the composition of the phytoplankton community evolves through time (Fasham et al., 1990). When one phytoplankton type becomes very abundant, specialized predators of that type also become more abundant and the integrated predator community becomes more efficient at grazing that phytoplankton type. Conversely, when a phytoplankton type is rare, its specialized predators will also be rare and the proportion of total grazing pressure on this phytoplankton type will be small.
From either perspective, the mathematical consequences are the same and so the parameterizations discussed in this study, including the KTW functional response, could be used to represent either ecological scenario. There is, however, an implicit assumption in use of switching to represent community-scale changes in the zooplankton. Namely, the functional response assumes that the composition of the zooplankton community changes instantaneously following the abundance of phytoplankton in the environment. This assumption is reasonable when used to represent microzooplankton predators, which typically have growth rates similar to their phytoplankton prey, however, it would be less appropriate when considering mesozooplankton predators that have significant generation times and seasonal-scale life cycles.
Functional responses that parameterize predator switching are an important tool in a modeler's toolbox, expanding the diversity of ecological phenomena that can be represented in food web models and providing a mechanism to promote diversity and counteract competitive exclusion. Work remains to be done to provide realistic constraints on the phenomenological parameters within the functional response, especially , which we have shown here has important dynamical consequences within plankton models that include KTW.

Appendix B. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jtbi.2023.111536.