INFERRING PARAMETERS OF PREY SWITCHING IN A PLANKTON 1 PREDATOR–2 PREY SYSTEM WITH A LINEAR PREFERENCE TRADEOFF

# Inferring Parameters of Prey Switching in a Plankton 1 Predator–2 Prey System With a Linear Preference Tradeoff

Sofia H. Piltz Corresponding author, Department of Applied Mathematics and Computer Science, Technical University of Denmark, Asmussens allé, Bygning 303B, 2800 Kongens Lyngby, Denmark and Department of Mathematics, University of Michigan, 2074 East Hall, Ann Arbor, Michigan, 48109-1043, USA (piltz@umich.edu).    Lauri Harhanen Department of Applied Mathematics and Computer Science, Technical University of Denmark, Asmussens allé, Bygning 303B, 2800 Kongens Lyngby, Denmark (lauri.harhanen@alumni.aalto.fi).    Mason A. Porter Department of Mathematics, University of California, Los Angeles, Los Angeles, California 90095, USA; Oxford Centre for Industrial and Applied Mathematics, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK; and CABDyN Complexity Centre, University of Oxford, Oxford, OX1 1HP, UK (mason@math.ucla.edu).    Philip K. Maini Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK and CABDyN Complexity Centre, University of Oxford, Oxford, OX1 1HP, UK. (maini@maths.ox.ac.uk).
###### Abstract

We construct two ordinary-differential-equation models of a predator feeding adaptively on two prey types, and we evaluate the models’ ability to fit data on freshwater plankton. We model the predator’s switch from one prey to the other in two different ways: (1) smooth switching using a hyperbolic tangent function; and (2) by incorporating a parameter that changes abruptly across the switching boundary as a system variable that is coupled to the population dynamics. We conduct linear stability analyses, use approximate Bayesian computation (ABC) combined with a population Monte Carlo (PMC) method to fit model parameters, and compare model predictions quantitatively to data for ciliate predators and their two algal prey groups collected from Lake Constance on the German–Swiss–Austrian border. We show that the two models fit the data well when the smooth transition is steep, supporting the simplifying assumption of a discontinuous prey switching behavior for this scenario. We thus conclude that prey switching is a possible mechanistic explanation for the observed ciliate–algae dynamics in Lake Constance in spring, but that these data cannot distinguish between the details of prey switching that are encoded in these different models.

##### Keywords:

prey switching, Lotka–Volterra interactions, linear stability analysis, planktonic ciliate–algae dynamics, smoothening out a discontinuous diet switch, parameter fitting using PMC–ABC

## 1 Introduction

Ciliates are eukaryotic single cells that propel using small protuberances (called cilia) that project from their cell body. They feed on small algae and are an important link between the bottom and higher levels of aquatic food webs [47]. In addition to seasonal variation, ciliates and their algal prey populations vary at shorter-than-seasonal temporal scales. During years when the spring bloom lasts for several weeks (corresponding to approximately 15–30 ciliate generations), algal and ciliate biomasses exhibit recurring patterns of growth followed by decline [47]. Ciliates have different modes of predatorial behavior, and they can be categorized roughly in terms of more-selective or less-selective feeding habits [53]. Some ciliate species can be construed as selective predators, because they hunt as “interception feeders” that scavenge on food particles and intercept them directly. In contrast, “filter-feeder” ciliates sieve suspended food particles and are an example of a less-selective ciliate species. A laboratory experiment on ciliate predator and phytoplankton prey species in Lake Constance reported prey preference and selective feeding in ciliates [33], and it has been suggested that predator–prey interactions between diverse predator and prey plankton communities are the driving force for the sub-seasonal temporal variability observed in ciliate–algal dynamics, especially during periods of the year in which environmental conditions are relatively stable [49].

In the present paper, we aim to obtain biological insight into the sub-seasonal oscillations in ciliate populations during spring in Lake Constance and more generally into the ecological concept of prey switching [34], in which predators express a preference, e.g., for more-abundant prey. To do this, we construct multiple modeling frameworks for a ciliate predator that adaptively changes its diet in response to changes in the abundances of its two prey.111“Prey switching” in a system of 1 predator and 1 prey refers to a situation in which predation is low at low prey densities but saturates quickly at a large value when the prey is abundant. In such a scenario, one can model the predator–prey interaction using a Holling type-III functional response [18, 12]. Using ciliate–algae interactions in Lake Constance as an example, we focus on adaptive feeding of a predator group between its two different types of prey to investigate both the characteristics of prey switching (specifically, whether it is best described with a steep or gradual switching function) and whether it is justified to use a reduced modeling framework (specifically, a piecewise-smooth dynamical system222Piecewise-smooth dynamical systems are a class of discontinuous systems that describe behavior using smooth dynamics of variables that alternate with abrupt events [9, 6].) as a simplifying approximation of a smooth system.

One can model prey switching with smooth dynamical systems by considering either density-dependent switching [2] or density-independent switching [37], or by using information on which prey type was last consumed [52, 51]. In contrast, a piecewise-smooth system arises when one assumes that a switch in a predator’s feeding behavior depends on prey abundances. For example, one can posit that a predator behaves as an optimal forager, as its choice to switch prey depends on which diet composition maximizes its rate of energy intake [42, 26, 30, 4, 28, 29]. Using a piecewise-smooth model, it was suggested recently that prey switching is a possible mechanistic explanation for the dynamics observed in ciliate and algae populations in Lake Constance [35].

In addition to their ecological applications, piecewise-smooth dynamical systems occur in a wide variety of applications [9], ranging from mechanical oscillators such as a rocking block (see, e.g, [16]) to relay-feedback systems (in which an electrically operated switch is used to control a process or an electromechanical system [10]). Other biological applications include gene regulatory networks, in which transcription factors either initiate or inhibit the production of proteins after some threshold concentration has been reached [13, 5], and conceptual climate models, in which an abrupt change in a piecewise-smooth system can represent a transition between different regimes in, e.g., large-scale ocean circulation [43] or the Earth’s reflectivity due to ice cover [1]. In a piecewise-smooth dynamical system, phase space is divided into two or more smooth regions by one or more switching manifolds that mark transitions between the regions. For prey switching, each region corresponds to one of the predator’s diet choices, so a model that describes the dynamics has a discontinuous right-hand side. Specifically, in this example, the system satisfies a different set of ordinary differential equations (ODEs) in different regions of phase space. In a piecewise-smooth framework, one assumes that a switch from one diet to another occurs instantaneously. Consequently, piecewise-smooth dynamical systems are also used to approximate nonlinear terms, such as sigmoidal or cubic functions, in models of systems that have sharp transitions between two or more states.

In both ecology and numerous other applications, it is important to conduct detailed investigations into different approaches for how to model sharp changes in governing dynamics. On one hand, it is unclear whether there exist “discontinuous predators” who switch their feeding strategy instantaneously, as assumed in a piecewise-smooth model for prey switching. On the other hand, we have not found evidence for any of the possible smooth transition functions that one can choose to model prey switching. By using data on ciliate and algal population dynamics, we illustrate an important example situation in which exploiting different modeling frameworks increases understanding of prey switching and allows one to gain biological insight into its profile. Consequently, we can conclude that it is justified to use a piecewise-smooth dynamical system, which has fewer parameters than the associated smooth models introduced in this paper, as a simplifying approximation of a smooth dynamical system. We also show that the piecewise-smooth model in [35] is both biologically and mathematically consistent as the limit of two smooth systems, which we construct by (1) using a hyperbolic tangent as a transition function from one diet choice to another and (2) incorporating a parameter that changes abruptly across the discontinuity in the model in [35] as a system variable with dynamics on a time scale comparable to that of the population dynamics of a predator and its two prey. In the second construction, we examine a system with one more dimension than the corresponding piecewise-smooth system.333Mathematically, one can derive different smooth approximations to the same piecewise-smooth system either by “smoothing out” a discontinuity of a piecewise-smooth system using a differentiable transition function of sigmoidal form [7] or by “regularizing” a piecewise-smooth dynamical system [27] into a singular perturbation problem [15, 24] that includes multiple time scales by “blowing up” the switching boundary [41, 44]. In this work, we do not consider regularizations that include multiple time scales. A subset of us were among the authors of previous work [36] on a multiple time-scale system describing the dynamics of one predator and two prey populations in the presence of rapid evolution of the predator’s diet choice. One can also include nonlinear terms when constructing a smooth dynamical system by smoothing out an instantaneous switch by using the method developed in [21, 22]. These nonlinear terms take into account small effects that are observable only during a switch and vanish for the corresponding piecewise-smooth system [23]. Comparing different smoothed out or regularized systems both with each other and with an associated piecewise-smooth system is crucial for understanding the correspondence and transition between a piecewise-smooth system and its smooth analogs. This type of relationship between piecewise-smooth and smooth systems that they approximate was investigated from a theoretical point of view in [17, 8].

The remainder of our paper is organized as follows. In Section 2, we present and briefly discuss the equations for the 1 predator–2 prey piecewise-smooth model from [35]. This piecewise-smooth dynamical system includes a tilted switching manifold that marks a transition between two smooth parts of phase space. Biologically, these parts represent the predator’s adaptive feeding behavior and its two different diet choices: on one side of the switching manifold, the predator’s diet consists solely of its preferred prey; on the other side, it consists solely of the alternative prey. We consider two possible regularizations of the model in Sections 2.1 and 2.2, and we derive analytical expressions and carry out linear stability analysis for the coexistence equilibrium (i.e., where all three species coexist at non-zero densities) for each of the two smooth models. We are interested in the coexistence steady states, because the data that we use includes coexistence of predators and multiple prey. In Section 3, we discuss and use these data on adaptively feeding plankton predators to fit model parameters and compare the biomass predictions of our two smooth models. We summarize and discuss similarities and differences in model behavior and model assumptions between the piecewise-smooth system (analyzed in [35]) and its two smooth analogs (analyzed in this paper) in Section 4, and we conclude our study in Section 5. We give additional details about our calculations and analysis in a trio of appendices.

## 2 Methods

We construct two smooth analogs of a piecewise-smooth dynamical system describing a predator population that can adjust the extent of its consumption of its preferred prey . When not consuming the preferred prey , the predator feeds on its alternative prey . Before introducing our two smooth models (see Sections 2.1 and 2.2), we present the model equations for the piecewise-smooth dynamical system that was developed in [35]. For each of the three models, we consider standard nonlinearities in the form of Lotka–Volterra predator–prey interactions. Although these nonlinearities are standard, it is convenient for us to use non-standard notation for the model coefficients that describe them. This notation allows us both to derive the switching condition introduced previously in [35] and to compare the two smooth models that we develop in the present paper to this piecewise-smooth system.

We assume that the predator switches to consume only an alternative prey when it maximizes its fitness by doing so. To describe this situation, Ref. [35] developed the following piecewise-smooth dynamical system:

 ˙x=⎡⎢⎣˙p1˙p2˙z⎤⎥⎦=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩f+=⎡⎢⎣(r1−β1z)p1r2p2(eq1β1p1−m)z⎤⎥⎦, if h=β1p1−aqβ2p2>0f−=⎡⎢⎣r1p1(r2−β2z)p2(eq2β2p2−m)z⎤⎥⎦, if h=β1p1−aqβ2p2<0, (2.1)

where and (with ) are the respective per capita growth rates of the preferred and alternative prey, and are the respective death rates of the preferred and alternative prey due to predation, is the predator per capita death rate per day, and is the predator conversion efficiency. The parameters and are nondimensional parameters that represent the predator’s respective desire to consume the preferred and alternative prey. Thus, the proportion of predation that goes into predator growth is given by for the preferred prey and by for the alternative prey. In other words, the model in Equation (2.1) preserves the idea of a predator’s benefit from consuming prey being in proportion to predation amount, as is the case in the standard Lotka–Volterra model. Specifically, we consider constant conversion efficiency as a fraction of for each prey, so the benefit from feeding on the preferred and alternative prey are represented by and , respectively. With , we emphasize the reduced benefit that the predator obtains from the alternative prey compared to its preferred prey. The difference in the benefit is also where the assumed tradeoff lies, as the alternative prey invests energy in building predator defense mechanisms and is thus a “less edible” prey compared to the preferred prey , which does not invest energy in predator defense mechanisms. Consequently, we assume that the growth rate of the preferred prey is larger than that of the alternative prey (i.e., ).

In our constructions (see Sections 2.1 and 2.2) of two smooth analogs of (2.1) and to facilitate our comparison between the piecewise-smooth and smooth systems, we take for simplicity. In doing so, we assume that the predator exhibits adaptive feeding behavior by adjusting its preference (rather than its attack rate) to the governing prey densities. The parameter corresponds mathematically to the slope of the tilted switching manifold between the two vector fields in (2.1). Biologically, is the slope of the assumed linear tradeoff in the predator’s preference for prey. The reader is referred to [35] for a biological justification of these model assumptions, analysis of the model (2.1), and inferred parameter values for data from Lake Constance.

In the present paper, we construct and carry out linear stability analysis of two novel (to our knowledge) smooth models for an adaptively feeding predator and its two prey. First, in Section 2.1, we formulate the model in (2.1) as a three-dimensional smooth dynamical system with hyperbolic tangent functions. Second, in Section 2.2, we construct a four-dimensional smooth analog of (2.1) by supposing that the desire to consume the preferred prey changes between and across the discontinuity in the piecewise-smooth system (2.1) as a system variable that changes along with the population dynamics.

### 2.1 Smooth model I

We construct a smooth analog (which we call our “smooth model I”) of (2.1) using a hyperbolic tangent as a transition function. This yields the following equations of motion:

 ˙p1 =r1p1−β1p1z(1+tanh(k(β1p1−aqβ2p2))2)≡f(p1,p2,z), ˙p2 =r2p2−β2p2z(1−tanh(k(β1p1−aqβ2p2))2)≡g(p1,p2,z), (2.2) ˙z =eq1β1p1z(1+tanh(k(β1p1−aqβ2p2))2)+eq2β2p2z(1−tanh% (k(β1p1−aqβ2p2))2)−mz≡l(p1,p2,z),

where determines the steepness of the transition function and thus of switches in the predator’s feeding behavior. Equation (2.2) incorporates Lotka–Volterra dynamics, and and (where ) can be construed, respectively, as one predator’s benefit from eating its preferred and alternative prey. In Section 3.1, we will infer values of that best fit data from a particular freshwater plankton system. The data were collected in Lake Constance between 1979 and 1999, were presented originally in [46, 47] and were subsequently analyzed further in several papers (e.g., [49, 48]). See Section 3 for a description of the data.

#### 2.1.1 Linear stability analysis of smooth model I

In this work, we are interested in a steady state of (2.2) with . We calculate

 f =p1(r1−z(1+tanh(k(p1−aqp2))2))=0 ⇒z(1+tanh(k(p1−aqp2))2)=r1 (2.3)

and

 g =p2(r2−z(1−tanh(k(p1−aqp2))2))=0 ⇒z(1−tanh(k(p1−aqp2))2)=r2, (2.4)

so

 l=(eq1p1−m)r1+(eq2p2−m)r2=0. (2.5)

We obtain the steady-state solution , where

 ~z =r1+r2 (eq1~p1−m)r1+(eq2~p2−m)r2 =0 (2.6) tanh(k(~p1−aq~p2)) =r1−r2r1+r2.

Taking the inverse hyperbolic tangent on both sides of the third equation in (2.6) results in linearly independent equations for and . We thereby obtain a unique steady state at

 ~p1 =aqm(r1+r2)+eq2r2arctanh(r1−r2r1+r2)ke(q1aqr1+q2r2), ~p2 =m(r1+r2)−eq1r1arctanh(r1−r2r1+r2)ke(q1aqr1+q2r2), (2.7) ~z =r1+r2.

All three population densities are positive at the steady state when

 k>k0=eq1r1arctanh(r1−r2r1+r2)m(r1+r2). (2.8)

We use the Routh–Hurwitz criterion [38, 19] to investigate the stability of the coexistence steady state (2.7).

###### Proposition 2.1.

If , then the steady state (2.7) is stable if and only if .

###### Proof.

See Appendix A. ∎

###### Proposition 2.2.

If , then there exists so that the steady state (2.7) is stable if and only if .

###### Proof.

See Appendix A. ∎

From these results, we see that when is large, which corresponds to a predator with a steep tradeoff in its prey preference (i.e., a small increase in specialization towards the preferred prey comes at a large cost in growth from feeding on the alternative prey), the coexistence steady state (2.7) is stable for all . When the prey switching is steep (i.e., when ), the coexistence steady state (2.7) is equal to the steady state of the piecewise-smooth system that lies on the switching manifold (i.e., it is a pseudoequilibrium), and it has a complex-conjugate pair of eigenvalues with negative real part when [35]. However, in contrast to the piecewise-smooth system, in which the coexistence steady state is repelling for shallow or flat prey preference tradeoffs (i.e., when ), the smooth system (2.2) has an interval of intermediate prey-switching slopes (see Equation (A.9) for the expression for ) for which the coexistence state is also stable for .

For population densities at the stable coexistence steady state, smooth model I predicts that the predator density is determined solely by the prey growth rates, so it is affected neither by the slope of the tradeoff nor by the steepness of the diet switch (see Equation (2.7)). For a nearly flat tradeoff (i.e., when is small), the stable coexistence steady-state solution for the preferred prey is at its minimum (when other parameter values are as given in the caption of Figure 1). However, if in addition to a mild tradeoff, the predator’s prey switching is also gradual (i.e., is very small), then has a high concentration at the stable equilibrium (see the left panel of Figure 1). This pattern of minimum and maximum values is reversed for the alternative prey: The steady-state concentration of the alternative prey has a large value when is small, except for very small , when the steady-state value of is small (see the right panel of Figure 1). Such behavior of the steady state (2.6) suggests that is a singular limit of smooth system I (2.2).

### 2.2 Smooth model II

To regularize the three-dimensional piecewise-smooth system (2.1) into a four-dimensional smooth system, we construct expressions for the temporal evolution of the predator’s trait to accompany the population dynamics of the predator and two prey. Biologically, we are assuming that the predator’s desire to consume its preferred prey undergoes either rapid evolution [11] or phenotypic plasticity [25], which are the two main forms of adaptivity in organisms. We will comment on these model assumptions in Section 4. We thereby turn the parameter that abruptly changes across the discontinuity in the piecewise-smooth model (2.1) (i.e., when and when ) into a system variable that changes in response to prey abundance on the same time scale as the population dynamics in a smooth dynamical system.

To ensure similarity with the piecewise-smooth model (2.1), we assume that no preference towards the preferred prey amounts to a feeding mode of consuming only the alternative prey (i.e., ) and that maximum preference towards the preferred prey amounts to a feeding mode of consuming only the preferred prey (i.e., ). We incorporate this assumption with a bounding function in the expression for the temporal evolution of the predator’s trait. From the condition for prey switching derived using optimal foraging theory [42] in [35], we impose that the rate of change of the mean trait value is proportional to . That is, we assume that the predator’s choice to switch prey depends on prey abundances and which diet composition maximizes its rate of energy intake [42]. For simplicity, we also assume exponential prey growth and linear functional response as in the piecewise-smooth system (2.1) [35]. We thereby obtain the following dynamical system for the population dynamics coupled with temporal evolution of the predator trait:

 dp1dt =g1(p1,p2,z,q)=r1p1−qp1z, dp2dt =g2(p1,p2,z,q)=r2p2−(1−q)p2z, (2.9) dzdt =g3(p1,p2,z,q)=eqp1z+e(1−q)q2p2z−mz, dqdt =f(p1,p2,q)=q(1−q)(p1−aqp2).

Similar to the piecewise-smooth system (2.1) and to smooth system I (2.2), the predator–prey interaction in (2.9) (which we call “smooth model II”) is of standard Lotka–Volterra type, so the benefit of consuming prey is proportional to the amount of predation. Consequently, the proportion of predation that goes into predator growth is given by for the preferred prey and by for the alternative prey, where is a parameter that represents conversion efficiency. For in the piecewise-smooth system (2.1), smooth model II (2.9) reduces to in (2.1) when and to in (2.1) when . Biologically, these two cases correspond, respectively, to the situations in which the predator’s diet is composed solely of the preferred prey and alternative prey. Note that the model in (2.9) does not include a time-scale difference, which was incorporated between demographic and trait (i.e., ) dynamics of the four-dimensional smooth system, and analyzed using singular perturbation theory in [36].

### 2.3 Linear stability analysis of smooth model II

The population densities of the predator and two prey at the coexistence steady state in smooth system II (2.9) are

 ~p1 =aqm(r1+r2)e(r1aq+r2q2), ~p2 =m(r1+r2)e(r1aq+r2q2), ~z =r1+r2, (2.10) ~q =r1r1+r2.

These are the same densities that occur both at the steady state of smooth system I (2.7) with steep prey switching (i.e., when ) and at the pseudoequilibrium point of the piecewise-smooth system (2.1) (for and ) that is located on the discontinuity boundary of the piecewise-smooth 1 predator–2 prey model [35].

We summarize the result of linear stability analysis of smooth system II (2.9) in the two following propositions.

###### Proposition 2.3.

If , then all eigenvalues of the steady state are purely imaginary.

###### Proof.

See Appendix B. ∎

###### Proposition 2.4.

If , then the steady state is linearly unstable.

###### Proof.

See Appendix B. ∎

Consequently, smooth system II (2.9) has an unstable coexistence steady state irrespective of whether a predator can be construed as selective with a steep preference tradeoff with respect to its preferred and alternate prey or as unselective with a mild tradeoff in its preference towards the two prey. Our results also imply that our smoothing of the piecewise-smooth system (2.1) by adding an extra dimension as in Equation (2.9) changes the stability of the coexistence equilibrium.

## 3 Results

To obtain insight into the steepness of the prey switching in the two smooth models that we constructed in Section 2, we consider data from Lake Constance (see [46, 48]) for ciliate predators and two different types of their algal prey groups. The Lake Constance data set consists of over 23,000 observations of abundances (expressed either as individuals or as cells per milliliter) and biomass (expressed as units of carbon per square meter) of various plankton species obtained at least once in a sample of a few milliliters to a liter of water between March 1979 and December 1999. We compare the abundances predicted by our two smooth models with data from years 1991 and 1998. (For a comparison between the piecewise-smooth model (2.1) and data, see [35].) During these two years, the spring bloom lasted for several weeks, and ciliate and algal biomasses exhibited recurring patterns of increases followed by declines [47]. We are interested in spring abundances, because previous studies have suggested that predator–prey feeding interactions are an important factor in explaining the ciliate–algae dynamics in that season [49] and that such interactions are more important than environmental conditions during spring [40].

Müller and Schlegel observed ciliates actively selecting against certain types of prey when offered a mixed diet of different types of their algal prey [33]. They suggested that adaptive feeding in ciliates occurs because different species benefit differently depending on the match between their feeding mode and the species that are abundant in the prey community. That is, ciliates select against their less-edible prey (e.g., a prey type that develops a hard silicate cover as a predator defense mechanism) when offered a mixed diet of both easily-digested and less-edible prey [33]. As a representative of an easily-digested prey group (i.e., preferred prey in the model equations), we consider data for Cryptomonas ovata, Cryptomonas marssonii, Cryptomonas reflexa, Cryptomonas erosa, Rhodomonas lens, and Rhodomonas minuta in the Lake Constance data set. For the less-edible prey (i.e., the alternative prey group ), we use data for small and medium-sized Chlamydomonas spp. and Stephanodiscus parvus. In addition to different prey groups, one can categorize ciliate-predators, which dominate the herbivorous zooplankton community in spring [47], roughly in terms of being more-selective or less-selective predators [53]. To represent differences in selectivity between different predator species, the unselective filter-feeder predator group consists of data for Rimostrombidum lacustris, and the selective interception-feeder predator group consists of data for Balanion planctonicum.

In this section, we use Lake Constance data on ciliate-predators and their algal prey to infer the steepness of the prey-switching function in smooth model I (2.2), the perturbation in the predator population from the coexistence steady state in Equation (2.10) of smooth model II (2.9), prey growth, predator death rates, and other parameters of our models. For our comparison between the Lake Constance data and the two smooth models in Sections 3.1 and 3.2, we first normalize both the biweekly data points and the model predictions for the predator density, , by their -norm (i.e., by Euclidean distance). We consider the time window from 1 March to 15 June, for which there are 31 data points for the selective predator and 19 data points for the unselective predator in 1991. In 1998, there are 15 data points for both the selective and unselective predator species between 1 March and 15 June. We fit parameters to data with approximate Bayesian computation (ABC) combined with a population Monte Carlo (PMC) method [3]. This combination allows us to study the results from the posterior parameter distribution rather than just a single value that gives the best fit as a result of an optimization method. Additionally, this approach allows us to code every step of the algorithm on our own. (We implement the algorithm in Matlab [45].) We can thereby examine possible sources of error in the fitting process more easily than if providing input and analyzing an output of an available program package for parameter estimation of ODE models. The posterior parameter distribution, which is an output of the fitting algorithm, is especially useful for assessing how well the piecewise-smooth model (2.1) approximates prey switching, which we represent with a hyperbolic tangent function in smooth model I (2.2) and by incorporating an additional system variable in smooth model II (2.9).

### 3.1 Comparison of smooth model I simulations and Lake Constance data

In this section, we fit the growth rates ( and , respectively) of the preferred and alternative prey, the predator mortality rate, , the slope of the prey-switching function, and the slope of the prey-preference tradeoff of smooth model I (2.2). We use as a bifurcation parameter. (See Propositions 2.1 and 2.2.) However, for simplicity (and similar to the study of the piecewise-smooth system in [35]), we assume that the nondimensional preference parameters are fixed (and we take and ). Thus, given our choice of the preference parameters and using as a bifurcation parameter, we investigate linear preference tradeoffs that all go through point but do so with different slopes.

Smooth model I (2.2) reproduces the peak abundances in the Lake Constance data and predicts an oscillatory pattern for both the selective and unselective predator populations during the springs of 1991 and 1998 (see Figures 2 and 3). Additionally, our parameter fitting suggests that adaptive feeding of the selective predator is best represented with a steep switching function. In particular, for 1998, we obtain more frequent gradual prey-switching functions for an unselective predator than for a selective one at the smallest tolerance level of the fitting algorithm. See the middle rows of Figures 2 and 3. Note that the coexistence steady state is unstable for the inferred parameter values that we use in the top-right panels in Figures 2 and 3. In these two figures, this steady state is thus unstable for and , respectively.

### 3.2 Comparison of smooth model II simulations and Lake Constance data

To compare simulations of the smooth model (2.9) to data, we fit the prey growth rates and , the predator mortality rate , and a perturbation in the predator population from the coexistence steady state in Equation (2.10) for (so that all four eigenvalues of the coexistence steady state are purely imaginary). We thus use as our initial value for the model simulations to infer values for that minimize the distance in Equation (C.1) between the data points and the model prediction for these points. Thus, a small perturbation suggests a gradual diet change and oscillating around its equilibrium value, whereas one can interpret a large perturbation from the equilibrium as rapid changes in the diet (and the dynamics of ).

Smooth model II (just like smooth model I) reproduces the peak predator densities, and it seems that smooth model II best fits the data when there is a large perturbation from the coexistence steady state. See Figures 4 and 5. For the year 1991, we predict that the selective predator group switches its diet less frequently than the unselective predator. Additionally, reaches its maximum and minimum values. In contrast, for the unselective predator, there is a change from increasing to decreasing in at some intermediate value. See the dynamics of at the bottom of the top panels of Figures 4 and 5. We predict that the switching of the selective predator occurs more often in year 1998 than in year 1991 (and in 1998, it occurs also more often than that of the unselective predator). See Figure 5.

To evaluate how well smooth model II (2.9) predicts prey abundance data (to which it was not fitted), we simulate it with parameter values that we obtain by fitting the model to the unselective predator in year 1991, and our model prediction for prey abundances suggests for an unselective predator in year 1991 that the preferred prey has smaller-amplitude oscillations than the alternative prey. As we show in Figure 6, this differs qualitatively from the data. We obtain the same result for smooth model I (2.2) (comparison not shown). Nevertheless, although we only use predator data to fit parameters, smooth model II (2.9) is also able to successfully capture some features of the prey data. As we illustrate in the left panel of Figure 6, such features include the periodicity of the peak densities of the preferred prey populations.

## 4 Discussion

From a biological perspective, using a smooth dynamical system allows us to relax the assumption of a “discontinuous” predator of the piecewise-smooth system (2.1). When the discontinuity is smoothed out using hyperbolic tangent functions, as in smooth model I (2.2), we can use data to determine the steepness of the transition in the predator’s feeding behavior for a particular predator type. Indeed, our parameter fitting to Lake Constance data suggests that one can model prey switching of either selective or unselective predator species with a steep hyperbolic tangent function. Additionally, our parameter fitting of smooth model I (2.2) predicts that the best fit to the data occurs in the parameter regime in which the coexistence steady state is unstable. Additionally, simulations of our smooth model II (2.9), which regularizes the abrupt change in the predator’s diet choice by considering a predator trait as a system variable, exhibit rapid predator-trait dynamics (i.e., the temporal evolution of the predator’s desire to consume the preferred prey ) suggesting that the best fit to data occurs when the change of diet is abrupt.

From a modeling perspective, the piecewise-smooth system (2.1) incorporates the effects of a predator’s adaptive change of diet in response to prey abundance, whereas smooth system II (2.9) (with an appropriate choice of parameters) explores rapid evolutionary change in a predator’s desire to consume its preferred prey [36]. Consequently, smooth system II models a different mechanism (namely, rapid evolution [11]) than the piecewise-smooth system (which models phenotypic plasticity [25]) for how rapid adaptation affects population dynamics [39, 54]. It has been suggested that one should be more likely to expect a stable equilibrium from models that account for phenotypic plasticity than in those that account for rapid evolution, because plastic genotypes respond faster than nonplastic genotypes to fluctuating environmental conditions [54]. Our modeling work is consistent with this hypothesis, as the piecewise-smooth system (2.1) converges to a steady state for a large region of phase space [35], but the same steady state is unstable — except for one point (), at which it is stable but nonhyperbolic — in smooth model II (2.9).

Smooth model I (2.2) necessitates the incorporation of a parameter that influences the system’s qualitative behavior, whereas smooth system II (2.9) has the same number of parameters as the piecewise-smooth system (2.1) but includes an additional system variable. It can thus be advantageous to study the piecewise-smooth system, especially if many species are included, because it allows one to avoid adding new parameters and/or variables. The hyperbolic tangent functions in (2.2) and the increased dimensionality of (2.9) both add complications to analytical calculations and parameter fitting. On the plus side, there are many more standard numerical techniques and standard theory to determine the stability of equilibria using linear stability analysis and study bifurcations for smooth dynamical systems than for piecewise-smooth ones. One needs to use more involved methods for theory and numerical computations for piecewise-smooth dynamical systems, and development of these techniques is an active area of research [9]. However, one can derive an analytical expression (specifically, ) for the flow at the discontinuity boundary of (2.1), and the available theory for piecewise-smooth dynamical systems identifies the bifurcation that takes place in (2.1) as crosses the value [35]. These results help facilitate understanding of the behavior of (2.1) and can be used when analyzing the ciliate–algae dynamics predicted by the piecewise-smooth model [35].

Both of our smooth models successfully reproduce the peak population densities and suggest that a parameter regime — when for smooth model I (2.2) and for a large perturbation from the steady state for smooth model II (2.9) — fits the data for ciliate predators in Lake Constance in the springs of 1991 and 1998. (Note that our initial distributions for these parameters in the fitting algorithm include both parameter regimes in which the coexistence steady state is stable and ones in which it is unstable.) Additionally, when using the parameters that we obtain from fitting smooth model II to data for the unselective predator in year 1991, we obtain agreement between our model’s prediction and the periodicity of the peak prey abundances. Both of our smooth models predict a higher frequency of peak densities than what we observe in the available data points. A high period in population oscillations is possible for small organisms, such as plankton, with short lifespans and large population densities. Making measurements more frequently than biweekly would be a good way to try to validate or refute the periodicity predicted by the smooth models. Additionally, using data comparison to choose between different competing models would be an effective way to help increase current understanding of the use of a piecewise-smooth model as a simplification of a steep transition in plankton-feeding behavior. Such comparisons would also be valuable more generally in numerous applications. In practice, one can carry out such a model comparison by implementing algorithms of model-based statistical inference methods (e.g., approximate Bayesian computation, as in the present paper) [50] or by using existing toolboxes for system identification (e.g., the ones implemented in Matlab [45]).

One can further investigate model predictions for the trait dynamics and compare them to results from controlled laboratory experiments with genetically diverse prey and/or predator populations in which one records the dynamics of the genetic diversity. Parameter fitting to Lake Constance data suggests that the best fit occurs in a parameter regime in which the predator trait dynamics oscillate abruptly between the maximum and minimum values. In a study of two plankton predators and their evolving algal prey, Hiltunen et al. [14] showed and discussed experimental evidence for periods of dominance of one predator followed by a rapid switch to a dominance by the other. In [14], the switch in predator dominance arose from interactions between changes in the predator population and changes in the frequency of the prey type that develops a predator defense mechanism against one of the two predators. Motivated by the above findings, it is also interesting to consider a model in which one incorporates a time-scale difference between demographic and predator-trait dynamics [36].

## 5 Conclusions

To increase biological insight into the experimentally-observed adaptive feeding behavior of unselective and selective ciliate predators on two different types of prey, we constructed two ordinary differential-equation models for prey switching. In one model (“smooth model I”), we represent the transition from one diet to another using the hyperbolic tangent function; in the other (“smooth model II”), we added a new system variable to describe the diet switch in the system (and we hence increased the system’s dimension by ). In constructing these models, we relaxed the simplifying assumption of a “discontinuous” predator feeding behavior in a piecewise-smooth dynamical system that was used previously to suggest prey switching as a putative mechanistic explanation for the observed dynamics [35]. Based on our results from fitting parameters of the two smooth systems to data on freshwater plankton, we conclude that the best fit to the data occurs when prey switching is rapid (and hence steep in the continuous models) and that the simplifying assumption of discontinuous predator feeding behavior appears to be justified.

Similar to earlier investigations, such as [20, 31], our study provides an illustrative example of both similarities and differences between a discontinuous system and smooth regularizations to it. When piecewise-smooth dynamical systems are used to simplify transitions in applications — such as approximating a cubic function in a membrane potential in models of spiking neurons [32], Hill functions in models of gene regulatory networks [13], changes in the Earth’s reflectivity due to ice melt in climate models [1], and more — understanding the extent to which the behavior of the corresponding smooth and piecewise-smooth systems agree is crucial for generating both accurate model simplifications and accurate predictions.

Finally, using the data that we currently possess, it is difficult to determine which of the three models (i.e., a piecewise-smooth model or the two smooth systems with an adaptive predator) provides a better mechanistic explanation for the observations of ciliate–algae dynamics in spring in Lake Constance. To enhance model selection, it would be very useful to collect data to improve analysis of the steepness of prey switching, the functional form of the preference tradeoff, and the periodicity of the population oscillations. Nevertheless, the construction of models using alternative modeling approaches, transforming between them, and comparing them to data can greatly increase understanding of the underlying mechanisms of biological systems.

## Acknowledgements

We thank John Hogan, Philip Maybank, Frank Schilder, and Frits Veerman for helpful discussions. We thank Ursula Gaedke for sending us the Lake Constance data, which were obtained as part of the Collaborative Programme SFB 248 funded by the German Science Foundation. SHP was supported by Osk. Huttunen Foundation (OHF) and Engineering and Physical Sciences Research Council (EPSRC) through the Oxford Life Sciences Interface Doctoral Training Centre and by People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement #609405 (COFUNDPostdocDTU). PKM would like to thank the Mathematical Biosciences Institute (MBI) at Ohio State University, for partially supporting this research. MBI receives its funding through the National Science Foundation grant DMS1440386.

## Appendix A Stability of the coexistence steady state in smooth model I

In this appendix, we prove Propositions 2.1 and 2.2.

###### Proof.

The Jacobian for Equation (2.2) is

 ⎡⎢ ⎢ ⎢ ⎢⎣r1−z2B(1+p1kC)zp1kaq2BC−p12Bzp2k2BCr2−z2C(1+p2kaqB)−p22Ceq1z2B(1+p1kC)−eq2p2zk2BCeq2z2C(1+p2kaqB)−eq1p1zkaq2BCeq1p12B+eq2p22C−m⎤⎥ ⎥ ⎥ ⎥⎦, (A.1)

where

 A =tanh(k(p1−aqp2)), B =1+A, (A.2) C =1−A.

At the coexistence steady state (2.7), the Jacobian (A.1) is given by

 1r1+r2⎡⎢⎣−2r1r2k~p12r1r2kaq~p1−r1~p12r1r2k~p2−2r1r2kaq~p2−r2~p2eq1r1(r1+r2)+2er1r2k(q1~p1−q2~p2)eq2r2(r1+r2)+2er1r2kaq(q2~p2−q1~p1)0⎤⎥⎦, (A.3)

which we henceforth denote by for the rest of the present appendix. The characteristic polynomial of (A.3) is

 |λI−J|=λ3+aλ2+bλ+c, (A.4)

where

 a=2k(p1+aqp2)r1r2r1+r2,b=1(r1+r2)2e(2kp21q1r21r2+p2q2r22(r1+2aqkp2r1+r2)+p1r1(−2kp2q2r1r2+q1(r21+r1r2−2aqkp2r22))),c=2ekp1p2r1r2(aqq1r1+q2r2)r1+r2. (A.5)

According to the Routh–Hurwitz criterion [38, 19], the coexistence steady state (2.7) is stable if and only if the coefficients in (A.5) satisfy , , and . The conditions and are satisfied because of the positivity of the system parameters. To study the third condition, we write as a polynomial in . We thereby obtain

 ab−c =2r1r2e(aqq1r1+q2r2)2(r1+r2)2(p1r1−aqp2r2)s(k),

where

 s(k)=s2k2+s1k+s0, (A.6)

with

 s0 =12e2q1q2r1r2log(r1r2)[2aqq1r1+2q2r2+(−aqq1r1+q2r2)log(r1r2)], s1 =em[(r1+r2)(a2qq21r21−q22r22)−r1r2(a2qq21r1+q22r2−3aqq1q2(r1+r2))log(r1r2)], s2 =4aqm2(aqq1−q2)r1r2(r1+r2).

Using the steady state (2.7), we see that

 p1r1−aqp2r2=aqkm(r21−r22)+e(aqq1+q2)r1r2arctanh(r1−r2r1+r2)ek(aqq1r1+q2r2)>0.

We have thus established that .

The value of at in Equation (2.8) is positive:

 (A.7)

For , the function is concave up with a positive derivative at , because

 s′(k0)=em(aqq1r1+q2r2)[(r1+r2)(aqq1r1−q2r2)+(3aqq1−q2)r1r2log(r1r2)]>0. (A.8)

Therefore, implies that for all . Finally, if , we see that is a downward-opening parabola. Because , it follows that is positive in the interval , where

 k1=−s1+√s21−4s2s02s2. (A.9)

## Appendix B Stability of the coexistence steady state in smooth model II

In this appendix, we prove Propositions 2.3 and 2.4. Both proofs use the characteristic polynomial of the Jacobian of smooth system II (2.9). This Jacobian is

 J =⎛⎜ ⎜ ⎜ ⎜⎝r1−~q~z0−~q~p1−~p1~z0r2−(1−~q)~z−(1−~q)~p2~p2~ze~q~ze(1−~q)q2~ze~q~p1+e(1−~q)q2~p2−me~p1~z−eq2~p2~z~q(1−~q)−aq~q(1−~q)0(1−2~q)(~p1−aq~p2)⎞⎟ ⎟ ⎟ ⎟⎠. (B.1)

At the coexistence steady state (2.10), the Jacobian (B.1) is given by

 J=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00−r1aqme(r1aq+r2q2)−aqm(r1+r2)2e(r1aq+r2q2)00−r2me(r1aq+r2q2)m(r1+r2)2e(r1aq+r2q2)er1er2q20e(aq−q2)m(r1+r2)2e(r1aq+r2q2)r1r2(r1+r2)2−aqr1r2(r1+r2)200⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (B.2)

whose eigenvalues are given by the roots of the characteristic polynomial

 det(λI−J)=λ4 +m(eq2r22+aqr1(er1+2r2))e(aqr1+q2r2)λ2 (B.3) +aqr1r2m2(aq−q2)(r1−r2)e(aqr1+q2r2)2λ+aqr1r2m2(r1+r2)e(aqr1+q2r2).
###### Proof.

(Proposition 2.3)

For , the term in Equation (B.3) vanishes. Substituting yields

 u2+m[2r1r2+e(r21+r22)]e(r1+r2)u+m2r1r2e. (B.4)

One can write the discriminant of (B.4) as

 D=m2e2(r1+r2)2⎛⎝(r21+r22)2(e−4r21r22(r21+r22)2)2+4r21r22(r21−r22)2(r21+r22)2⎞⎠. (B.5)

Note that is always positive, so the two roots ( and ) of (B.4) are both real. Furthermore, because the polynomial (B.4) is increasing and positive at the intersection of (B.4) with the vertical axis, the roots of the polynomial (B.4) are both negative. Consequently, the four eigenvalues (with ) consist of two complex-conjugate pairs with real part: and . We thus see that all eigenvalues are purely imaginary.

###### Proof.

(Proposition 2.4)

First, we prove by contradiction that there is at least one eigenvalue with a non-zero real part. Assume that all four eigenvalues are purely imaginary. One can then write the characteristic polynomial (B.3) as

 χ(λ)=4∏j=1(λ−ıyj),yj∈R. (B.6)

Expanding , we see that the coefficient is

 ı(y1y2y3+y1y2y4+y1y3y4+y2y3y4),

which is purely imaginary. However, (B.3) has a real coefficient for that is non-zero for . Therefore, there exists at least one eigenvalue with a non-zero real part.

To complete the proof, we show that there are two eigenvalues whose real parts have opposite signs. We denote the roots of (B.3) by (with ), and we order the roots so that the real part of is non-zero. Because , at least one of , , or must have a real part whose sign is opposite to that of . Consequently, the steady state is unstable. ∎

## Appendix C Parameter fitting

We study the parameter-fitting problem through Bayesian inference. In contrast to least-squares fitting, this allows one to study the results from the posterior parameter distribution rather than just a single value that gives the best fit as a result of an optimization method. Because our modeling results in a posterior associated to the data, we fit parameters to data with approximate Bayesian computation (ABC) combined with a population Monte Carlo (PMC) method. (See p. 987 of [3].)

### c.1 Smooth model I

Let denote the solution of (2.2) with initial values , and let denote the available measurement data on the predator population. The data were measured at time instances , so — without measurement errors — the data would be for some unknown, true parameter values. We account for the presence of measurement errors by incorporating normally-distributed noise into the results of our model simulations. Specifically, for given parameter values , the model prediction is described element-wise as , where is the maximum predator density in a model trajectory. We compute a model trajectory by simulating the model for about 400 days and discarding the first about 60 days (corresponding to the two winter months January and February) as a transient. We then align the peak abundances in the data and in the model trajectory obtained by simulating the model with the given parameter values. One can construe this procedure as introducing a phase shift in the model results before calculating the distance between it and the data 444Such a procedure results in several candidate parameter sets that need to be rejected (e.g., because they predict a steady state); this increases computation time. An alternative to using a phase shift is to fit the initial values simultaneously with the model parameters. In such an approach, one can compare the distances between the periodic orbits that result from the model to those in the data. However, it is not clear how one should choose a reasonable time window for fitting the initial values and whether such a modification would yield more effective parameter fitting than with the current approach.. Because we do not know the variance of measurement errors in advance, we incorpo