{'fraction_non_alphanumeric': 0.0401861442618029, 'fraction_numerical': 0.018269158597913383, 'mean_word_length': 4.332266068998639, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A predator-2 prey fast-slow dynamical system for rapid predator evolution

We consider adaptive change of diet of a predator population that switches its feeding between two prey populations. We develop a novel 1 fast--3 slow dynamical system to describe the dynamics of the three populations amidst continuous but rapid evolution of the predator's diet choice. The two extremes at which the predator's diet is composed solely of one prey correspond to two branches of the three-branch critical manifold of the fast--slow system. By calculating the points at which there is a fast transition between these two feeding choices (i.e., branches of the critical manifold), we prove that the system has a two-parameter family of periodic orbits for sufficiently large separation of the time scales between the evolutionary and ecological dynamics. Using numerical simulations, we show that these periodic orbits exist, and that their phase difference and oscillation patterns persist, when ecological and evolutionary interactions occur on comparable time scales. Our model also exhibits periodic orbits that agree qualitatively with oscillation patterns observed in experimental studies of the coupling between rapid evolution and ecological interactions.


Introduction
Organisms can adapt to changing environmental conditions-such as prey availability, predation risk, or temperature-by changing their behavior. For example, in prey switching, a predator changes its diet or habitat in response to prey abundances. This is an example of phenotypic plasticity [37], in which the same genotype can express different phenotypes in different environments. However, adaptivity can also be expressed as an evolutionary change in traits (i.e., properties that affect how well an individual performs as an organism [49]) via genomic changes of a predator and/or prey [22]. If such evolution occurs on a time scale of about 1000 generations and can be observed in laboratory conditions, it is construed as a 'rapid' evolutionary change of traits [22]. Rapid evolutionary changes have been observed in a wide variety of organisms, ranging from mammals [53] to bacteria [6], and both in predators (e.g., in traits that involve resource consumption [26] or the ability to counteract prey defense mechanisms [28]) and in prey (e.g., in traits that involve predator avoidance [35,72]).Understanding the dependencies between rapid evolution and ecological interactions is fundamental for making accurate predictions on a population's ability to adapt to, and persist under, changing environmental conditions [15,57]. For example, rapid evolutionary change of traits has been observed in a plankton predator-prey system [21,72], which is a good example system for studying the coupling between rapid evolution and predator-prey interaction due to its short generation times and the tractability of genetic studies of it [35].
Adaptive change of feeding behavior can be incorporated into a dynamical system of predator-prey interaction in multiple ways. For example, one can represent prey switching with a Holling type-III functional response [33], consider the densities of the different prey as system variables [2,56], or use information on which prey type was last consumed [66,67]. Such formulations lead to smooth dynamical systems, but one can also model a predator that switches prey using a piecewise-smooth dynamical system [7,14] in which continuous temporal evolution of predator and prey populations alternates with abrupt events that correspond to points at which the predator changes its diet or habitat [41,45].
Rapid evolution. Several theoretical and empirical investigations have considered the effect of rapid evolutionary change of traits on predator-prey dynamics (see [22] for a review), including examples in which ecological and evolutionary dynamics have been assumed to occur on (1) comparable time scales or (2) disparate time scales. An example of (1) is the occurrence of out-of-phase cycles between small zooplankton (i.e., the predator) and genetically variable clonal lines of algae (i.e., the prey) populations observed in the experiments in [21,72], which were reproduced using a mathematical model with contemporary evolutionary and ecological dynamics [36]. The model in [36] suggests that the cycles emerge from prey evolution, especially when there is a small (energy) cost associated with the prey defense mechanism. Examples of (2) include situations in which evolutionary change occurs on either a slower [39] or on a faster [9] time scale than ecological interactions. Consequences that ecological dynamics can have on trait evolution have also been studied using the mathematical framework of adaptive dynamics [23,73], where evolution is assumed to occur on a slower time scale than ecological interactions. In the present paper, we aim to give insight into how the evolution of traits arises in population dynamics, and we hence concentrate on studying the limit in which trait evolution occurs on a much faster time scale than predator-prey interactions. When the time scales can be separated, one can use the framework of fast-slow dynamical systems [44] to introduce and exploit a time-scale separation between evolutionary and ecological dynamics to reduce the dimensionality of the system of equations that describe the evolutionary and ecological interactions.
When evolutionary change is faster than ecological interactions, the fast-slow dynamical system introduced in [9] can preserve the qualitative properties of dynamics in a predator-evolving prey model in which ecology and evolution occur on the same time scale [36]. In addition, similar to a model with only one time scale [36], a fast-slow dynamical system with rapid prey evolution reproduces experimentallyobserved out-of-phase predator-prey oscillations [72]. However, such oscillations are not present in the analogous model without rapid evolution [9]. By exploiting the time-scale separation between fast evolution and slow ecology, the general theory of either an evolving prey or an evolving predator [9] has been extended to cover the case in which both predator and prey evolve [8]. There now exist general conditions for determining which type of cyclic dynamics are possible in a system of coevolving prey and predator [8]. Such dynamics involve cycles that exhibit (i) counterclockwise or clockwise orientation in the predator-prey phase plane, (ii) a half-phase difference between the predator and prey oscillations, and (iii) "cryptic" cycles in which the predator population cycles while the prey population is approximately constant. Interestingly, a situation in which both predator and prey evolve can generate clockwise cycles that have also been identified in empirical data sets from systems such as phage-cholera, mink-muskrat, and gyrfalcon-ptarmigan [10]. This contrasts with traditional Lotka-Volterra predatorprey cycles which have a counterclockwise orientation in the phase plane with a quarter-phase difference between the predator and prey oscillations [51].
Our approach. In the present paper, we use a similar approach to [8,9] and develop a novel (to our knowledge) fast-slow dynamical system for a predator switching between two groups of prey species. In contrast to [8,9], we make a simplifying assumption of unlimited prey growth (e.g., because of favorable environmental conditions) and use Lotka-Volterra functional responses between predator and prey. As a result, we can prove that there exists a family of periodic orbits in a system of one predator and two different prey species. As we discuss in Section 7, some of our orbits agree qualitatively with patterns observed in both laboratory experiments and observations. In addition to the potential utility of mathematical modeling (and using fast-slow systems) for understanding the coupling between ecology and evolution [22], our motivation for constructing our model comes from our earlier work that suggests adaptive prey switching of a predator as a possible mechanistic explanation for patterns observed in data on freshwater plankton [55]. In the model in [55], the switch in the predator's feeding behavior is discontinuous. In the present paper, we relax this assumption and consider a rapid but continuous change in the predator's feeding choice. 1 There exists theory both for regularizing a given piecewise-smooth system to create a fast-slow system [60,62] and for approximating a fast-slow system using a piecewise-smooth system that preserves-both qualitatively and quantitatively-key characteristics (such as singularities and bifurcations) of the original fast-slow system [12,13]. However, we instead construct our model from a biological perspective using a concept from quantitative genetics. A clear understanding of the trajectories of solutions of the fast-slow system makes it possible to compare the fast-slow system as an "ecologically-obtained" regularization (which we construct using fitness-gradient dynamics [3,46]) of the piecewise-smooth system with a regularization obtained from a mathematical viewpoint (e.g., using a method based on a blow-up technique [60,62]). We aim to shed light on the differences and similarities that phenotypic plasticity and rapid evolution, as two different mechanisms, can generate in a model for the population dynamics of an adaptively feeding predator and its two prey. In Section 7, we compare and contrast the dynamics exhibited by the earlier piecewise-smooth model [55] and the fast-slow system obtained as an ecological regularization of it.
Putting aside our motivating applications in ecology, we note in passing that classical examples of fast-slow dynamical systems include the Van der Pol [65] and Fitzhugh-Nagumo equations [19,52]. The former was used originally to describe the dynamics of an electrical circuit with an amplifying valve (and has subsequently been used for numerous other applications), and the latter is a simplified version of the Hodgkin-Huxley nerve-axon model [32] from neuroscience [17]. There are several other applications of such multiple time-scale systems-including pattern formation [24], opening and closing of plant leaves [20], ocean circulation [40], critical transitions in climate change [43], and more. We also note that multiple time-scale systems can be studied using several different techniques from singular perturbation theory [31,44]. Examples of such techniques include matched asymptotic expansions [5,38] and geometric singular perturbation theory [34,44]. Because we are interested in constructing periodic orbits and understanding their bifurcations in the fast-slow dynamical system for prey switching (see Section 2), we use the latter to analyze our model.
Geometric singular perturbation theory. We aim to gain insight into how the evolution of traits that occurs on a comparable time scale to that of ecological interactions arises in population dynamics by studying the limit in which trait evolution occurs on a much faster time scale than that of the predator-prey interactions. Consequently, our goal in choosing a geometric approach is to create solution trajectories for a parameter ε > 0 by concatenating segments of curves that are determined by either the fast reduced dynamics or the slow reduced dynamics when ε = 0. Following the notation in [44], a fast-slow dynamical system with m fast variables and n slow variables (and time as the only independent variable) is expressed as where f : R m × R n × R → R m , g : R m × R n × R → R n , and ε (with 0 < ε 1) is the ratio of the two time scales. We rescale the slow time t by ε and obtain an equivalent system that evolves on the fast time scale τ = t/ε. We thus write We can take the (singular) limit ε → 0 in Eq. (1.1), which describes the dynamics evolving on the slow time scale t, to obtain the reduced slow vector field We call (1.3) the slow reduced problem. Similarly, we can take the limit ε → 0 in Eq. (1.2), which describes the dynamics evolving on the fast time scale τ , to obtain the reduced fast vector field which is called the fast reduced problem. The fast and slow reduced problems are connected through the critical manifold C 0 = {(x, y) ∈ R m × R n : f (x, y, 0) = 0}, a sufficiently smooth submanifold of R m × R n .
The critical manifold C 0 determines the equilibrium points of the fast reduced problem (1.4), and the differential-algebraic slow reduced problem (1.3) determines a (slow) dynamical system on C 0 . By using geometric singular perturbation theory, we study the critical manifold and the slow (1.3) and fast (1.4) reduced problems to obtain information about the behavior of the full system (1.1). This approach builds on the work of Fenichel [18], which guarantees that, under some general conditions, several geometric objects (e.g., the critical manifold C 0 ) defined in the reduced slow and fast problems persist for sufficiently small ε > 0 as similar geometric objects in the full system (1.1). Hence, for example, (singular) orbits that are constructed by concatenating orbit pieces from the slow reduced problem (1.3) and the fast reduced problem (1.4) persist for sufficiently small ε > 0 in the sense that such a singular orbit is a good approximation of a "true" orbit of the full system (1.1).
For an introduction to geometric singular perturbation theory and its concepts, see [29]. For a comprehensive overview of singular perturbation theory, see [44].
Outline of our paper. The rest of our paper is organized as follows. In Section 2, we formulate a 1 fast-3 slow dynamical system for a predator feeding on two prey populations in the presence of rapid predator evolution. The three slow variables of the system correspond to the populations of the predator and the two prey. We model a predator trait that represents the predator's feeding choice as the fast variable of the system. This model construction allows us to use geometric singular perturbation theory to gain insight into the effects on population dynamics of an evolutionary change of a predator trait that occurs on a time scale that is comparable to that of the predator-prey interaction. In Section 3, we derive expressions for the critical manifold and the slow and fast subsystems. We then use these results in Section 4 to explicitly construct periodic orbits that are exhibited by the 1 predator-2 prey fast-slow system. We obtain these expressions for the periodic orbits by studying the singular limit in which the ratio ε of the fast to the slow time scale goes to 0. In Section 5, we highlight some ecologically relevant qualitative aspects of the constructed periodic orbits. We then use numerical continuation in Section 6 to investigate how the periodic orbits persist for ε > 0 as we perturb the system. Finally, we discuss the findings of our study in Section 7.

The fast-slow 1 predator-prey model
We begin our formulation of a 1 predator-2 prey fast-slow system in Section 2.1 by constructing an equation for the temporal evolution of a predator population (z) that adaptively changes its diet between two prey populations (p 1 and p 2 ). Our fast-slow model is based on four principal assumptions. We assume that the organisms (1) have a large population size, (2) live in a well-mixed environment, and (3) can be aggregated into groups of similar species. Consequently, we can represent the predator-prey interaction as a low-dimensional system of ordinary differential equations. We also assume that (4) it is possible for evolutionary changes of traits to occur on a time scale that is comparable to that of ecological interactions. In Section 2.2, we define an expression for the temporal evolution of a predator trait (q) that represents the predator's desire to consume each prey. In our model, we assume that the predator trait undergoes rapid evolution on a faster time scale than that of the population dynamics. In Section 6, we use this assumption to gain insight into interactions between ecological and evolutionary dynamics when these occur on a comparable time scale.

Ecological dynamics
We assume that the predator's desire q to consume prey is bounded between its smallest and largest feasible values (q S and q L , respectively). For simplicity, we consider exponential prey growth and a linear functional response between the predator growth and prey abundance [51]. We thereby obtain the following system of differential equations for the population dynamics of the 1 predator-2 prey fast-slow system: where r 1 and r 2 (with r 1 , r 2 > 0) are the respective per capita growth rates of the prey p 1 and p 2 , the parameters β 1 and β 2 are the respective death rates of the prey p 1 and p 2 due to predation, e > 0 is the proportion of predation that goes into predator growth, q 2 ∈ [0, 1] is the nondimensional parameter that represents the extent of preference towards prey p 2 , and m > 0 is the predator's per capita death rate. One can also interpret the parameter q 2 as a factor that scales the benefit that the predator obtains from feeding on prey p 2 . For simplicity, we let β 1 = β 2 (which we can take to be equal to 1 by rescaling) to omit β 1 and β 2 in our calculations in Section 3. In doing so, we assume that the predator exhibits adaptive diet choice by adjusting its feeding choice (i.e., whether the predator is feeding on prey p 1 or on prey p 2 ) rather than its attack rate based on the prey densities. We also require that the extreme when q is at its minimum (i.e., q = q S ) corresponds to the case in which the predator is feeding solely on prey p 2 . Similarly, we require that the extreme when q is at its maximum (i.e., q = q L ) corresponds to the case in which the predator feeds solely on prey p 1 . We thereby assume that q is bounded between q L and q S . Without loss of generality, we choose q L = 0 and q S = 1. These assumptions simplify the system (2.1), which represents the fast-slow 1 predator-2 prey population dynamics, to dp 1 dt = r 1 p 1 − qp 1 z ,

Evolutionary dynamics
We assume that the adaptive change in the predator's trait q follows fitness-gradient dynamics. In other words, we assume that the rate of change of the mean trait value is proportional to the fitness gradient of an individual with this mean trait value [3]. Fitness-gradient dynamics was used for defining trait dynamics of rapid predator evolution in [8][9][10], wherein fast-slow dynamical systems were proposed as a general framework for gaining insight into evolutionary and ecological dynamics that occur on a comparable time scale. In its original form in [3], the fitness F of an individual is assumed to be frequency-dependent. That is, F = F (q * , q), where q * is the trait value of an individual and q is the mean trait value of the population. In the present paper, we determine fitness as the net per capita growth rate of the predator population. The rate of change of the mean population trait value is then governed by the fitness-gradient dynamics as follows: Following ecological considerations, if the adaptation occurs by genetic change, then the rate constant that describes (2.3) is the additive genetic variance [46] (i.e., genetic variance due to genes whose alleles contribute to the trait value additively). For simplicity, we assume that the fitness of an individual depends only on the mean trait value of the population and not on the distribution of the individual's trait. That is, F (q, q * ) = F (q). In addition, we describe the additive genetic variance of the predator's desire to consume each prey type by a bounding function that limits the predator trait between its smallest (q S = 0) and largest (q L = 1) feasible values. Furthermore, by assuming that the predator trait evolves on a faster time scale than the population dynamics (where the separation of time scales is given by ε), we see that the temporal evolution of the predator trait takes the following form: where V is a nondimensional constant and is part of the additive genetic variance term q(1 − q)V . When using fitness-gradient dynamics, we model evolutionary dynamics at the phenotypic level without incorporating detailed information about genotypic processes (e.g., principles of Mendelian inheritance [22]). Using this simplified approach, we can incorporate an equation for the predator trait directly into the system of 1 predator-2 prey dynamics and obtain an analytically tractable differential-equation model for coupled ecological and evolutionary dynamics. Because of its simplifying assumptions on genotypic processes and the laws of inheritance, fitness-gradient dynamics gives an incomplete understanding of interactions between ecological and evolutionary dynamics [22]. Nevertheless, there is evidence that fitness-gradient dynamics can still be an appropriate approximation for modeling evolutionary dynamics even when its simplifying assumptions do not hold [3,25].

Coupled ecological and evolutionary dynamics
By combining the ecological dynamics in (2.2) with the evolutionary dynamics in (2.4), we obtain the following fast-slow 1 predator-2 prey system with predator evolution: When q = 1 in (2.5), the predator feeds only on prey p 1 . Likewise, when q = 0, the vector field of the fast-slow system (2.5) corresponds to a situation in which the predator's diet is composed solely of p 2 . Consequently, there is exponential growth in the population of the prey type that is not being preyed upon.

Analytical setup
In this section, we introduce concepts from geometric singular perturbation theory to aid in the analysis of the 1 fast-3 slow model in (2.5).

Rescaling of the system (2.5)
To keep our analysis as clear as possible, we rescale Eq. (2.4) to maximally reduce the number of parameters. Using the rescaling we obtainṗ where r and m are free parameters. Without loss of generality, we can assume that 0 < r < 1 (i.e., that r 1 > r 2 ).

Linearization around the nontrivial equilibrium of (3.2)
The system (3.2) has a unique nontrivial equilibrium point at (p 1 , p 2 , z, q) = 1, 1, 1 + r, 1 1+r . The local behavior near the equilibrium is characterized by the eigenvalues of the linearization of (3.2) around it. These eigenvalues (λ) obey the characteristic equation For all m > 0 and 0 < r < 1, equation (3.3) has four purely imaginary solutions, so the nontrivial equilibrium is of center-center type. For an equilibrium of this type, local analysis alone is in general extremely complicated and intricate (see, e.g., Ch. 7.5 of [27]). Moreover, any local periodic solution stays close to the value q = 1 1+r , whereas our goal is to study the effect of a prey-switching mechanism in which q varies between-and should come near-the extremal values 0 and 1. We therefore abandon the standard linearization approach and turn to geometric singular perturbation theory to obtain farfrom-equilibrium periodic orbits (see Section 4).

Analysis of system (3.2) as a fast-slow system
In this section, we specify the fast and slow subsystems of the full system (3.2) and compute the critical manifold C 0 . We use the symbol " · " to denote derivatives with respect to the slow time t and and use the symbol " " to denote derivatives with respect to the fast time τ .

Slow reduced system
We obtain the slow subsystem of the full system (3.2) by considering the singular limit ε → 0. In this limit, we obtain the slow reduced systeṁ As we described in Eq. (1.3), this is a differential-algebraic system.

Fast reduced system
We scale the slow time t in Eq. (3.2) by ε and reformulate its dynamics in terms of the fast time τ = t/ε to obtain Taking the limit ε → 0 yields p 1 = 0, p 2 = 0, and z = 0. The reduced fast system is one-dimensional, and it determines the fast dynamics of q through We can solve Eq. (3.6) explicitly for q(τ ) because p 1 and p 2 are constant. For 0 < q < 1, we obtain which is gauged such that q h (0) = 1 2 .

Critical manifold
The critical manifold C 0 is defined by the algebraic part of the slow reduced system (3.4). It is given by where We see that the critical manifold can be written as the union of a trio of three-dimensional hyperplanes. Observe that the hyperplane M 0 (3.9) is an invariant manifold for the slow reduced system (3.4). Indeed, any initial condition with q(0) = 0, when evolved according to (3.4), yields q(t) = 0 for all t. This allows us to study the flow of the slow reduced system (3.4) on the hyperplane M 0 through the dynamical systemṗ The dynamics of p 1 decouples from the variables p 2 and z, and the prey p 1 exhibits exponential growth. The predator z and the prey p 2 form a Lotka-Volterra predator-prey system around the nontrivial equilibrium (p 2 , z) = (1, r). We can thus introduce a conserved quantity on M 0 ; it is given by 3.3.5 Slow flow on the hyperplane M 1 (3.10) By the same reasoning as in Section 3.3.4, the hyperplane M 1 is an invariant manifold for the slow reduced system (3.4). The flow of (3.4) on M 1 is given bẏ These dynamics are very similar to those on the hyperplane M 0 (3.12), but the roles of p 1 and p 2 are reversed. Now, the dynamics of p 2 decouples from the variables p 1 and z, and the prey p 2 exhibits exponential growth. The predator z now forms a Lotka-Volterra predator-prey system with the prey p 1 around the nontrivial equilibrium (p 1 , z) = (1, 1). The associated conserved quantity on M 1 is given by

Construction of periodic orbits
In this section, we use the setup from Section 3 to provide a geometric analysis of Eq. (3.2) in terms of its slow (3.4) and fast (3.5) subsystems. We then use results from geometric singular perturbation theory to combine these two descriptions to construct a family of (singular) periodic solutions to (3.2) for sufficiently small ε. We also indicate how one can construct these solutions explicitly using the analytical results in Section 3. In Section 6, we show the numerical simulations for the constructed periodic solutions for specific values of ε.

Geometric analysis of the fast reduced system
As we explained in Section 1, the critical manifold consists of equilibrium points of the reduced fast system (3.6). Geometrically speaking, the flow defined by (3.6) connects a point (p 1 , p 2 , z, 0) ∈ M 0 with a point (p 1 , p 2 , z, 1) ∈ M 1 by a heteroclinic connection, which is given explicitly by Eq. (3.7). The sign of p 1 − p 2 determines the direction of this heteroclinic connection. Therefore, the hyperplane M sw divides the four-dimensional phase space into two parts (see Figure 1): on the side in which p 1 > p 2 , the heteroclinic connection q h (3.7) is directed from M 0 to M 1 (i.e., "upwards"); on the other side of M sw , in which p 1 < p 2 , the direction of the heteroclinic connection is reversed, going from M 1 to M 0 (i.e., "downwards"). The hyperplane M sw acts as a switching plane: when this plane is crossed, the direction of the heteroclinic flow that connects M 0 and M 1 is reversed. Figure 1: Flow of the fast reduced system (3.6). For visual clarity, we depict the three slow model dimensions (p 1 , p 2 , z) in two dimensions, which are spanned by (p 1 − p 2 , z); the vertical axis indicates the fast variable q. We indicate the hyperplane M 0 in green and the hyperplane M 1 in blue. In this visualization, the switching hyperplane M sw is spanned by the z-axis and q-axis. On the right side of M sw , the fast flow (indicated by double arrows) is directed upwards (i.e., from M 0 to M 1 ); on the left side of M sw , the direction of the fast flow is reversed (i.e., it is directed downwards).

Geometric analysis of the slow reduced system
We study the flow of the slow reduced system (3.4) on M 0 , as given by (3.12), from a geometric point of view. The phase space of the flow on M 0 is three-dimensional, and it is given in terms of the coordinates (p 1 , p 2 , z). Projected onto the (p 2 , z)-plane, the system (3.12) reduces to a classical Lotka-Volterra system with the conserved quantity H 0 (3.13). We know [51] that every orbit of this Lotka-Volterra system is closed and is determined uniquely by its value of H 0 . Because the p 1 -dynamics are decoupled from the (p 2 , z)-dynamics, every H 0 level set in the (p 2 , z)-plane extends to a cylindrical level set in the full (p 1 , p 2 , z) phase space. Because H 0 is also a conserved quantity for the full system (3.12), these cylindrical level sets of H 0 are invariant under the flow of (3.12). Therefore, we can characterize the dynamics of (3.12) by describing the three-dimensional phase space as a concentric family of cylindrical level sets of H 0 (see Figure 2). For the flow of the slow reduced system (3.4) on M 1 , as given by (3.14), a geometric perspective yields an equivalent construction, with the roles of p 1 and p 2 reversed. In this case, one can characterize the dynamics of (3.14) on the same (p 1 , p 2 , z) phase space through another concentric family of cylindrical level sets, which are determined by the conserved quantity H 1 (3.15) (see Figure 3).

Combining the fast and slow reduced dynamics
We seek to use the geometric insight in the fast and slow reduced limits of the full system (3.2) (see Consider Figure 4). On this side of the switching plane M sw , the heteroclinic connection q h (3.7) takes us "up" to the corresponding point We now take the viewpoint of the slow reduced limit, and we use the point A 1 = (p A 1 , p A 2 , z A , 1) as an initial condition for the slow flow on M 1 . Let Φ t 1 denote the slow flow on M 1 . We let the slow flow on M 1 act for some time T 1 , so the point A 1 ∈ M 1 flows to the point For notational brevity, we write Φ T1 1 p A 1 = p B 1 , and we use similar notation for p A 2 and z A , writing B 1 = (p B 1 , p B 2 , z B , 1). We choose T 1 so that the p 1 -coordinate is now smaller than the p 2 -coordinate (i.e., p B 1 < p B 2 ). Switching back to the viewpoint of the fast reduced limit, we see (because p B 1 < p B 2 ) that the slow flow on M 1 has brought us to the other side of the switching plane M sw . We can thus use the heteroclinic connection q h to travel "down" from B 1 ∈ M 1 to the corresponding point B 0 = (p B 1 , p B 2 , z B , 0) ∈ M 0 . Back on M 0 , we again take the viewpoint of the slow reduced limit, and we take the point B 0 = (p B 1 , p B 2 , z B , 0) as an initial condition for the slow flow Φ t 0 on M 0 . We let the slow flow on M 0 act for some time T 0 . Because our goal is to construct a periodic (and hence closed) orbit, we want to choose (p A 1 , p A 2 , z A ) and the times T 0 and T 1 so that the slow flow on M 0 takes B 0 back to the starting point A 0 .
We give a schematic overview of the above construction in the following diagram:

Existence of singular periodic orbits for ε > 0
We successfully constructed a singular orbit in Section 4.1 by concatenating several orbit sections that we obtained by studying the reduced slow (3.12), (3.14) and fast (3.6) limits of the full system (3.2).  Using geometric singular perturbation theory, which builds on work by Fenichel [18], one can show that when ε is nonzero but small, orbits constructed in this singular setting "persist". In other words, for ε > 0 but sufficiently small, the full dynamical system (3.2) includes periodic orbits that are very close to their singular counterparts. The following theorem makes this observation more precise. Because these orbits in the full system (3.2) still exhibit a sharp distinction between "slow" and "fast" segments, it is conventional to call them "singular".
Theorem 4.1 (Persistence of singular periodic orbits). Let ε > 0 be sufficiently small, and let γ 0 (t) be a singular closed orbit as constructed in Section 4.1. Then there exists a periodic orbit γ ε (t) in the full system (3.2) that is O(ε) close to γ 0 (t).
Proof. One can prove the existence of the periodic orbit γ ε by following the procedure in Soto-Treviño [59]. However, as the reduced fast system (3.6) is one-dimensional, the slow manifolds M 0 and M 1 lose normal hyberbolicity at the intersection with M sw (for which p 1 − p 2 = 0), where the system exhibits a slow passage through a transcritical bifurcation [42]. Near M 0 ∩ M sw , we use the standard blowup transformation t = √ εt, p 1 − p 2 = √ εp, q = √ εq to obtain, up to O( √ ε), the equations d dt (p 1 + p 2 ) = 0 = dz dt and dp dt It is clear that there is no exchange of stability [47,48]. There is one canard, which is maximal and is given byq = 0. Furthermore, using the results in [11], it follows that 0 <q(t) <q(−a) for all −a <t < a with a ∈ O(1). In other words,q(t) stays close toq = 0 for O(1) time int. The analysis at M 1 ∩ M sw is similar. Combining these observations with the analysis in [59] yields the existence of γ ε in an O(ε) neighborhood of γ 0 . For more details on the proof of the construction, see [59]; for more details on the local analysis near p 1 − p 2 = 0, see [11].

Existence conditions and solution families
The goal of this section is to establish (algebraic) conditions for the existence of the closed singular orbit that we described in Section 4.1.3. The periodicity of the orbit in the fast q-coordinate is satisfied by construction through the heteroclinic connections between M 0 and M 1 . Therefore, we only need conditions on the initial values of the slow coordinates (p 1 , p 2 , z) and on the slow-evolution times T 0 and T 1 . We need to choose these five unknowns in a way that ensures periodicity in the three slow coordinates. Therefore, we generically expect to obtain a (possibly empty) two-parameter family of periodic orbits 3 . See Figure 6 for a visualization of the singular orbit that depicts the relevant quantities that we use in the following analysis. The expression (3.13) for the conserved quantity H 0 (p 2 , z) of the reduced slow flow on M 0 provides a relation between the slow coordinates of the "take-off" point A 0 and the slow coordinates of the "touch-down" point B 0 . We obtain the relation from Likewise, on M 1 , we use the expression (3.15) for the conserved quantity H 1 (p 1 , z) to obtain a relation between the slow coordinates of the touch-down point A 1 and the take-off point B 1 . We obtain the relation from We now consider the explicit form of the slow reduced dynamics on M 1 . On M 1 , we flow the point A 1 to the point B 1 using the flow (3.14) for a time T 1 . The flow of the p 2 -coordinate is decoupled from the variables p 1 and z, and it is linear, so we solve for its dynamics directly to obtain The other two slow coordinates (p 1 , z) interact through Lotka-Volterra dynamics. We integrate Eq. (3.14) forṗ 1 to obtain where we can obtain the expression z(p 1 ) by invoking the conserved quantity H 1 and inverting the relation . That is, we obtain z(p 1 ) by solving the equation One can treat the slow segment on M 0 analogously. On M 0 , we flow the point B 0 back to the point A 0 using the flow (3.12) for a time T 0 . The flow of the p 1 -coordinate is decoupled from the variables p 2 and z, and it is also linear, so we can solve for its dynamics directly to obtain On M 0 , the slow coordinates (p 2 , z) interact through Lotka-Volterra dynamics. We integrate the equation forṗ 2 using (3.12) to obtain where we obtain the expression z(p 2 ) using the conserved quantity H 0 by inverting the relation H 0 (p 2 , z) = H 0 (p A 2 , z A ). That is, we obtain z(p 2 ) by solving the equation We can use the above results on the slow flow on M 0 and M 1 to eliminate T 0 and T 1 . Combining Eq. (4.4) with (4.5) yields 1 and combining Eq. (4.7) with (4.8) yields Together with (4.2) and (4.3), we have now obtained four equations for six unknowns, which consist of the slow components of A 0,1 and B 0,1 . The relations (4.6) and (4.9), which need to be inverted to obtain the integrands of (4.5) and (4.8), are not bijective. Therefore, the specific forms of z(p 1 ) (4.5) and z(p 2 ) (4.8) depend on the characteristics of the underlying slow-orbit segment. To obtain computable expressions for the coordinate values of (p A 1 , p A 2 , z A ) and (p B 1 , p B 2 , z B ), these underlying slow-orbit segments need to be characterized in more detail. We demonstrate this procedure in Appendix A. We have numerically evaluated the explicit expressions that we thereby obtain for the coordinate values of (p A 1 , p A 2 , z A ) and (p B 1 , p B 2 , z B ) for several choices of the model parameters r and m. In the Online Supplement, we show the results of these numerical evaluations, together with visualizations of the associated singular periodic orbits.

Ecologically relevant qualitative aspects of the constructed periodic orbits
In this section, we discuss several qualitative aspects of periodic solutions that are ecologically relevantincluding synchronization between predator and prey and/or between two prey species, clockwise cycles, and counterclockwise cycles. These types of behavior occur (1) in the data collected from microscopic aquatic organisms in the field [63,64] and under laboratory conditions [4,30,71] and (2) in experimental studies of coevolution in phage-bacteria systems [50,69]. Importantly, such behavior also arises in the family of periodic orbits that we constructed in Section 4.

Prey-prey synchronization
We have access to field data on microscopic aquatic organisms [63,64], and these data exhibit both antiphase and in-phase oscillations between the two different prey types. Although there are several experimental studies of coevolution in microscopic aquatic organisms (and in phage-bacteria systems that, as another kind of exploiter-resource system, translate to a model of predator-prey interaction), we have not yet encountered empirical observations of evolutionary and demographic dynamics in a system of one predator and two different prey species. Note, however, that in a bacterium-phage system studied in [68], a bacterial (i.e., prey) subpopulation that replicates slowly and is phenotypically resistant to the phage (i.e., predator) was suggested as a possible mechanistic explanation for the observed dynamics.
This agrees with our model assumption of the growth rate of the alternative prey being less than that of the preferred prey (i.e., 0 < r < 1 in the rescaled system in Eq. (3.2)). When two populations oscillate in (anti)phase, local extrema of two species occur at exactly the same instances in time. Because p 1 increases monotonically on M 0 and p 2 increases monotonically on M 1 , it follows that p 1 does not have an extremum during the slow dynamics on M 0 and that p 2 does not have an  (p 1 , z); on the right, we show the phase plane spanned by (p 2 , z). In each phase plane, we use a solid curve to indicate the orbit segment in which the displayed model variables are interacting through Lotka-Volterra dynamics. We indicate the associated level curves of (left) H 1 (p 1 , z) (3.15) and (right) H 0 (p 2 , z) (3.13) in grey. We use a dashed curve to indicate the remaining orbit segment, in which the displayed prey variable grows exponentially. We use black dots to indicate the jump points A and B, at which the periodic orbit "jumps" from M 0 to M 1 , and vice versa. The arrows on the orbit segments give the direction of time.
extremum during the slow dynamics on M 1 . Therefore, for p 1 and p 2 to oscillate in (anti)phase, we need their local, aligned extrema to occur at the jump points-i.e., where q jumps from 0 to 1, or vice versa. Because p 1 increases from jump point B towards jump point A, it follows that p 1 has a local maximum at A and a local minimum at B. By the same reasoning, p 2 must have a local minimum at A and a local maximum at B. Therefore, we can conclude that the only type of prey-prey synchronization that occurs in the singular periodic solutions that we have constructed occurs when the two prey species oscillate in antiphase. Therefore, based on the above considerations, the dual phase-plane picture associated with prey-prey synchronization must be as depicted Figure 7. As is clear from this figure, the z-coordinates of both jump points A and B must lie above the z = 1 nullcline. That is, z A > 1 and z B > 1. Our numerical calculations show that there exists a two-parameter family of periodic orbits in which both prey species oscillate in antiphase for a range of the model parameters (r, m); see Figure 8 for a visualization of such a solution family. See the Online Supplement for visualizations of the associated singular periodic orbits and for our numerically obtained values of (p A 1 , p A 2 , z A ) and (p B 1 , p B 2 , z B ).

Predator-prey synchronization
Because our model includes two prey species, predator-prey synchronization can potentially arise either through synchronization between the predator and prey 1 or though synchronization between the predator and prey 2. Predator and prey densities that oscillate almost exactly out of phase with each other have been observed in experimental studies on the effects of rapid prey evolution on ecological dynamics [4,72].
In this section, we examine the conditions under which there is a jump from one slow manifold to the other, and we thereby show that our model exhibits (anti)phase oscillations in which the predator is in phase with a prey at one jump point and out of phase with the same prey at the other jump point.
Synchronization between all three species. We observe that the predator cannot oscillate either in phase or in antiphase with both prey species simultaneously, because this would imply that the aligned local extrema of the prey species would be of the same type (i.e., both maxima or both minima). In other words, the two prey species would exhibit in-phase oscillations, and we showed in Section 5.1.1 that this cannot occur. However, it is still possible for the predator to oscillate in phase with one prey and in antiphase with the other prey. Suppose that the two prey oscillate in antiphase, as described in Section 5.1.1. If the predator is in phase with one prey and in antiphase with the other, then the predator density has local extrema located at the jump points. The nature of these extrema is dictated by the "jump conditions" on p 1,2 at the jump points A and B in the following manner. Suppose that the predator density z has a local maximum at A. From the slow flow on M 0 (3.12) and the slow flow on M 1 (3.14), we infer that this implies that (p A 2 − 1)mz A > 0 and (p A 1 − 1)mz A < 0. However, this violates the jump condition that p A 1 > p A 2 (see Section 4.1.3). Therefore, z cannot have a local maximum at A. By an analogous argument, we see that z cannot have a local maximum at B. Therefore, if the predator density has local extrema at both jump points A and B, then both of these extrema are local minima. We conclude that the only way in which the predator is synchronized with both prey species is if the predator is alternating between (1) being in phase with prey 1 and in antiphase with prey 2 and (2) being in antiphase with prey 1 and in phase with prey 2. For an example of such predator-prey-prey synchronization, see Figure 9.
Synchronization between predator and one prey. Suppose that the predator z oscillates in (anti)phase with one prey only, which we assume is p 1 without loss of generality. From the slow dynamics of p 1 on M 1 (3.14) and on M 0 (3.12), it is clear that during the slow dynamics, the local extrema of z and p 1 are unable to align. On M 0 , the prey p 1 changes monotonically; on M 1 , the prey p 1 and predator z are related through Lotka-Volterra dynamics, which forbids alignment of local extrema of the participating predator and prey. We thus conclude that the local extrema of the predator z and the prey species p 1 with which it oscillates in (anti)phase must occur at the jump points A and B. From the fact that p 1 increases monotonically during the slow dynamics on M 0 , we conclude that p 1 must be a maximum if it has a local extremum at A. This, in turn, implies that the derivative of p 1 on M 1 at the jump point A must be negative. Therefore, (1 − z A )p A 1 < 0 (see Eq. (3.14)), so z A > 1. For the other prey species p 2 , we see using (3.14) that the derivative of p 2 on M 1 at A is (trivially) positive and that the derivative of p 2 on M 0 at A is given by (r −z A )p A 2 . However, we have already concluded that z A > 1, and because we have also assumed that 0 < r < 1 (see Section 3.1), the derivative of p 2 on M 0 at A must be negative. Therefore, if p 1 has a local extremum at A (which must be a maximum), then p 2 has a In the panels of this figure, we show the projections of these six-tuples onto different coordinate planes; see also Figure 16.   local minimum at A. This violates the assumption that the predator z oscillates in (anti)phase with one prey only. We therefore conclude that the predator z cannot oscillate in (anti)phase with prey 1 only. Any synchronization between the predator and prey 1 necessarily implies synchronization between the predator and prey 2. We discussed this predator-prey-prey synchronization in the previous paragraph.
It is worthwhile to note that the simultaneous synchronization of both prey is a direct consequence of the asymmetry between the dynamics of p 1 and p 2 , manifested through the parameter 0 < r < 1 in (3.14) and (3.12). Now suppose that the predator z oscillates in (anti)phase with one prey only, and suppose that that prey is p 2 . By the same arguments as above, we can conclude that if p 2 has local extrema at A and B, then p 2 must have a local minimum at A and a local maximum at B. Therefore, from (3.14) and (3.12), it follows that r < z A and r < z B . Because we assumed that z does not oscillate in (anti)phase with prey p 1 , the derivative of p 1 must be positive on both sides of jump point A, and it must also be positive on both sides of jump point B. This implies that z A < 1 and z B < 1. We can conclude that it is possible for the predator z to synchronize with one prey only, and this prey must be prey 2. This situation occurs if and only if r < z A < 1 and r < z B < 1. However, from the previous paragraph, we know that the local extrema of the predator z at the jump points A and B can only be minima. Therefore, the predator z and prey p 2 cannot oscillate in phase or in antiphase. The only synchronization possible between z and p 2 is of an "alternating" type: if z and p 2 are synchronized, then they are in phase at jump point A and in antiphase at jump point B. For an example of this type of synchronization, see Figure 10.
Based on the above analysis, we see that the only phase-locking mode available for the periodic orbits that we have constructed in this paper is a "hybrid" phase in which the predator is in phase with one prey at one jump point and out of phase with the same prey at the other jump point. In summary, the above analysis on synchronization types yields the following insights.
1. Alignment of local extrema of the model species (p 1 , p 2 , z) can occur only at a jump point.
2. If the predator z has a local extremum at a jump point, then this extremum must be a local minimum.
3. When z A > r and z B > r, the predator z and prey p 2 are in phase at A and in antiphase at B.

4.
When z A > 1 and z B > 1, the predator z and prey p 1 are in antiphase at A and in phase at B.
5. Therefore, the only two types of synchronization between the predator z and prey p 1,2 are the ones in Figures 9 and 10.

Clockwise and counterclockwise cycles
In this section, we discuss the ordering of the peak abundances of the predator and prey populations in cycles exhibited by the model in Eq. (3.2). In particular, we describe both situations in which a peak in the predator abundance precedes that in the prey population (so that the cycles have a "clockwise" orientation when depicted on a predator-prey phase plane), and those in which a peak in the prey abundance is followed by that in the predator population (so that the flow travels "counterclockwise" on the phase plane). The nomenclature "clockwise" and "counterclockwise" stems from the orientation of the flow in a classical Lotka-Volterra system, which describes one prey species and one predator species. In the traditional phase-plane depiction of the Lotka-Volterra system, the prey is placed on the horizontal axis, and the predator is placed on the vertical axis. In this case, the Lotka-Volterra flow has a counterclockwise orientation. In the solution time series, this dynamical behavior is characterized by the fact that a peak in the prey population is relatively close to a peak in the predator population. Moreover, in Lotka-Volterra dynamics, the prey peaks first and the predator peaks shortly thereafter.
The prediction of counterclockwise cycles due to density-dependent predator-prey interactions in the Lotka-Volterra model is supported by empirical observations collected from hare and lynx populations [16]. As is also the case for several other traditional predator-prey models, the Lotka-Volterra system assumes that the behavior and characteristics of the organisms remain fixed on the time scale of ecological interactions. As we discussed in Section 1 and will discuss further in Section 5.1, rapid evolution alters the population dynamics and, in particular, it can generate cycles in which the peak in the predator abundance follows the peak in the prey population with a phase lag that is larger than a quarter of a period [4].
In contrast to the counterclockwise cycles, clockwise cycles are characterized by a negative phase lag between the peak abundances and a reversed ordering of the predator and prey maxima. Recently, Cortez and Weitz [10] analyzed ecological data sets collected from various different predator-prey systems and identified regions in them that have a clockwise orientation. In [10], a peak in the predator population is considered to precede a peak in the prey abundance if the time between a predator peak and the following prey peak is less than the distance between the predator peak and the preceding prey peak. Theoretical modeling work suggests that evolutionary changes on a time scale comparable to that of the ecological interactions and occurring in both predator and prey offer a possible mechanistic explanation for the reversed ordering of the peak abundances [10]. In other words, a small population of a prey type that has invested in predator defense mechanisms and a large population of a predator that is ineffective in counteracting the prey's defense can yield low predator abundance and high prey abundance because of the effective prey defense. Consequently, selection favors predators that are effective in counteracting prey defense so that the prey population starts to decrease. Simultaneously, the predator population remains low because of the high cost of counteracting prey defense. However, due to low predator population, there is room for the prey population with low predator defense to increase. The predator population then increases because of a high abundance of more vulnerable prey [10]. In this reversed situation, the predator peaks first, and the prey peaks shortly thereafter.
In the predator-prey-prey system that we study in the present paper, the answer to the question of whether clockwise cycles occur seems to be straightforward. Using the orbit visualization in Figure  7, we immediately see that both the interaction between the predator and prey 1 and the interaction between the predator and prey 2 occur in a counterclockwise fashion. Such orientations are inherent to the construction of the (singular) periodic orbits in question (see also Figure 5). However, from such phase portraits, one is unable to draw any conclusion about the difference in time between predator and prey extrema. In particular, using a phase-space perspective, one cannot readily deduce whether a prey peak is shortly followed by a predator peak (the time-series hallmark of a counterclockwise cycle), or vice versa.
To obtain more insight in the relative time difference between predator and prey peaks, we use the analysis of Section 5.1. The slow dynamics on M 0 [see Eq. (3.12)] and M 1 [see Eq. (3.14)] show that on either slow manifold, one prey species increases monotonically and the other prey species interacts with the predator through Lotka-Volterra dynamics. In Lotka-Volterra dynamics, a peak in prey density always precedes a predator peak. Therefore, the situation in which a predator peak is shortly followed by a prey peak cannot occur during the slow dynamics on either M 0 or M 1 . However, as we saw in Section 5.1, local prey maxima can also occur at the jump points A and B. (The predator can only have local minima at the jump points; see Subsection 5.1.1.) In Figure 11, we show an example of a singular periodic orbit in which the prey density p 1 has a local maximum at A, where q jumps from 0 to 1. The predator peak occurs shortly before this instance in time. In Figure 12, we show another  Figure 11: Rescaled abundances of the preferred prey p 1 (solid blue curve), alternative prey p 2 (dashed blue curve), the predator z (dotted red curve), and predator trait q (dot-dashed black lines) as a function of the rescaled time t (on the horizontal axis) for a singular orbit exhibiting "clockwise" behavior (i.e., the peak in the predator density occurs just before the peak in prey 1) in the system in Eq. singular periodic orbit, in which the predator peak occurs almost exactly in between the peak of prey 1 and the peak of prey 2. We can therefore conclude that the dynamics of our model (2.5) admits (singular) periodic solutions whose time series exhibit an ordering of predator and prey peaks that can be interpreted as "clockwise" in the sense that a prey peak is shortly preceded by a predator peak. In Section 6, we will show that the localization in time of these local maxima persists when we increase the value of the asymptotically small parameter ε.

Numerical continuation of the singular periodic orbits
In this section, we use direct numerical simulation of the model system (3.2) to illustrate the theoretical results from Section 4. The goal of this section is to highlight the role of the small parameter ε. We demonstrate that we can numerically find the singular periodic orbits that we constructed in Section 4, which we proved to exist for ε 'sufficiently small' (see Theorem 4.1). We also demonstrate that these periodic orbits persist as ε is increased to larger values (even ones for which we no longer have a theoretical guarantee that such an orbit exists). In other words, we perform a numerical continuation in ε starting from a singular periodic orbit in which ε = 0.
To initiate the numerical-continuation procedure, we use the explicit analytical solution of the singular periodic orbit, which is characterized by the slow coordinates of the jump points A and B. For a specific choice of (r, m), we find the values of (p A 1 , p A 2 , z A ) and (p B 1 , p B 2 , z B ) that satisfy the existence conditions (4.2), (4.3), (4.10), and (4.11). We consider a singular periodic orbit, for which a peak in the predator population lies between the peaks in the two prey populations (for a singular orbit of this kind, see Figure 12), and we use the parameter values (r, m) = (0.5, 0.4). After several numerical simulations, we obtain (p 1 , p 2 , z, q) ≈ (1.18, 0.87, 1.50, 0.99) and use these values as an initial condition, simulate the system (3.2) with ε = 0.025, and find a numerical solution that is nearby the corresponding singular orbit (see panel (a) of Figure 13). We carry out the continuation of this solution for increasing values of ε as follows. We simulate the rescaled system (3.2) (for 50 time units when 0.025 ≤ ε ≤ 0.5 and for 30 time units when 0.5 ≤ ε ≤ 1) and use the final value of each simulation as an initial value for the next simulation in 3 sequences of 10 simulations with ε linearly spaced between 0.025 ≤ ε ≤ 0.2 (for the simulation with ε = 0.2, see panel (b) of Figure 13), 0.2 ≤ ε ≤ 0.5 (for the simulation with ε = 0.5, see panel (c) of Figure 13), and 0.5 ≤ ε ≤ 1 (for the simulation with ε = 1, see panel (d) of Figure 13).
One can clearly observed that, as the value of ε increases, the transition in q between its minimum and maximum values is increasingly gradual. Moreover, in our numerical simulations, we see that the theoretical minimal value of q (i.e., q = 0) is not attained by the solutions in panels (b), (c), and (d) of Figure 13. From numerical continuation, we see that choosing the small parameter ε = 0.025 allows us to numerically find periodic orbits that are O(ε) close to the singular limit. That is, for the parameter choice (r, m) = (0.5, 0.4), the value ε = 0.025 can be construed as "sufficiently small" in Theorem 4.1. In addition, we see that certain quantitative features of the singular periodic orbit persist as ε is increased. These include the antiphase oscillation between the two prey species and the placement of the predator peak in between the prey peaks. Thus, the analytical results obtained through geometric singular perturbation theory, which are established for sufficiently small values of ε (see Theorem 4.1), are seen to be meaningful for 'unreasonably large' values of ε.

Conclusions and discussion
Summary and time scales. We have modeled adaptive feeding behavior of a predator that switches between two prey types. In our model (see Eq. (2.5)), we assume that the predator gradually changes its diet from one prey to another depending on the prey densities and that the predator feeds only on one prey type at a time at the extremes. The change of diet is continuous, but it is fast compared to the time scale of population dynamics, so we have introduced a time-scale difference between the dynamics of a predator trait (which represents the predator's desire to consume each prey) and the population dynamics. The resulting 1 fast-3 slow dynamical system exhibits periodic orbits that we can construct analytically when the parameter ε that represents the separation of time scales is equal to 0, and we have proven that these orbits persist for small but nonzero values of ε. We have also demonstrated, using numerical computations, that these periodic solutions persist for (non-small) values of ε from ε = 0.025 up to ε = 1.
Beyond piecewise-smooth formulations: various ways to smoothen a jump. Part of our motivation to construct a fast-slow dynamical system for a predator adaptively switching between two prey comes from the desire to relax the assumption of a "discontinuous" predator that was made in earlier work [55]. The discontinuity in this previous model, which successfully reproduces the periodicity in the ratio between the two prey groups exhibited in data collected from freshwater plankton [55], comes from the assumption that a predator chooses a diet that maximizes its growth [61]. Because it is not clear whether there exist predators that switch their feeding strategy instantaneously, our future work includes two other ways in which we "smooth out" the 1 predator-2 prey piecewise-smooth system in [54]. One can regularize a piecewise-smooth dynamical system into a singular perturbation problem by "blowing up" the discontinuity boundary. This method was developed originally in [60] and later surveyed in [62]. In the present paper, we combined these ideas with the ecological concept of fitness-gradient dynamics [3,46] to relate this discontinuity smoothening to a biological phenomenon. Fitness-gradient dynamics was used previously to represent trait dynamics of a predator-prey interaction in a fast-slow system for studying rapid evolution and ecological dynamics [9]. Note, however, that the fast-slow system in [9] is not obtained as a result of regularizing a given piecewise-smooth system. Instead, the rate of change of a predator (or prey) trait is assumed to be governed by fitness-gradient dynamics [3,46] and to evolve at a faster time scale than that of the predator-prey interaction [9] (as in the present paper).
Obtaining analytical results using time-scale separation. An important feature of our model (3.2) is the presence of the parameter ε, which introduces a time-scale separation between the (fast) model component q and the (slow) model components (p 1 , p 2 , z). The inclusion of this time-scale separation enables us to not only include a biological mechanism (i.e., natural selection) directly into the model, but also to rigorously prove the existence of periodic orbits for sufficiently small values of this time-scale separation parameter ε. Using numerical methods, we have demonstrated that these "sufficiently small" values of ε are within numerical reach, as we are able to find values for ε for which direct numerical simulation of the model (3.2) yields periodic orbits that are O(ε) close to their singular counterparts. By using numerical continuation, we have also shown that these orbits persist for increasing values of ε. Therefore, the method that we outlined in this paper can be used not only to study periodic solutions to (3.2) in the presence of a time-scale separation, but also (using numerical continuation of the singular periodic solution) to investigate the existence and behavior of periodic solutions to (3.2) when the time scale of the rate of change of the trait q is comparable to that of the predator-prey interaction. The periodic solutions that we constructed in this paper are far-from-equilibrium solutions, so the model variables donot stay close to the system's nontrivial equilibrium. As we mentioned in Section 3.2, this equilibrium is of center-center type. The use of local analysis around this equilibrium to study periodic solutions is in general extremely complicated and intricate, and the nature of this type of analysis excludes the study of far-from-equilibrium solutions (and, in particular, it excludes ones in which the trait variable q switches between 0 and 1. It is worth noting that the method that we employed in this paper to construct (singular) periodic orbits is not confined to orbits with two slow segments-one on each slow hyperplane M 0,1 , as in Figure 4. Using the same methods, our analysis can be extended to study periodic orbits with two slow segments on each slow hyperplane by concatenating them using four fast transitions. This would lead to an extension of the family of possible periodic orbits. This larger class of periodic orbits, which have various qualitative features, could be fitted to experimental data.
Rapid evolution versus phenotypic plasticity. The two mechanisms of adaptivity-i.e., phenotypic plasticity and rapid evolution-both cause rapid adaptation and affect population dynamics [58,70].
Although it is not clear precisely how these different mechanisms affect population dynamics, it has been suggested that models that account for phenotypic plasticity exhibit a stable equilibrium more often than models that account for rapid evolution [70]. It has also been suggested that this can arise from a faster response time of plastic genotypes than that of nonplastic genotypes to fluctuating environmental conditions [70]. Indeed, our model (2.5)-which describes the population dynamics of a predator and its two prey in the presence of rapid evolutionary change in a predator trait-does not contain stable steady states. In contrast, the model in [55]-which considers an adaptive change of diet in response to prey abundance-exhibits convergence to a steady state for a wide parameter range.
Cryptic and out-of-phase cycles. Empirical evidence suggests that rapid evolution is a possible mechanistic explanation for cyclic dynamics that differ from those in traditional predator-prey systems.
For an evolving prey, such dynamics include (1) large-amplitude cycles in a predator population while a prey population remains nearly constant [6,71], (2) predator and prey oscillating almost exactly out of phase [4,72], and (3) oscillations in which a peak in a prey population follows that in a predator population [4]. The types of oscillations (1) and (2) also arise in predator-prey models with rapid predator evolution [9]. It has also been demonstrated that rapid predator evolution as a response to an evolving prey can generate cyclic dynamics both in experiments [30] and in models [8] of coevolution. Our model (2.5) represents an evolving predator that feeds on two different types of prey and exhibits different types of periodic orbits-including ones in which the predator and prey populations oscillate out of phase, total prey density remains approximately constant, and a peak in the prey population follows that in the predator populations.
Future work and comparison with experiments and field observations. To identify orbits that exist in an ecologically reasonable parameter range, our ongoing work includes comparing our model simulations with data collected from freshwater plankton in the field [63,64] for which our principal model assumptions (i.e., large population size, short generation times, and well-mixed environment) hold. However, similar models can be formulated for any other organisms that satisfy these assumptions, including the microorganisms used in the laboratory experiments in [30]. The insight into rapid evolution gained from studying a tractable plankton system can be used as an example for understanding rapid evolution of larger organisms and their abilities to adapt to changing environmental conditions (such as climate change or species introductions). Moreover, one of the major applications of the understanding of coevolution in microorganisms is resistance to antimicrobial drugs. By fitting parameters of prey growth rates and predator mortality to data, we expect to be able to distinguish parameter regimes to determine which members of periodic-orbit families best describe the data, and one can thereby gain insights into a system of one evolving predator that feeds on two different types of prey. In particular, we expect such comparisons between models and data to help determine ecological trade-offs and their possible influence on rapid predator evolution. Empirical evidence from a study of coevolving predator and prey suggests that a predator pays a low fitness cost (or no cost at all) for counteracting anti-predatory prey evolution [30]. However, in addition to the unknown mechanism of the predator response, the trade-off(s) that constrain rapid predator evolution remain unknown if one only looks at data without doing any modeling.