Frequency‐dependent evolution in a predator–prey system

We present a discrete‐time predator–prey model in which the prey population is assumed to evolve in response to a toxicant. We incorporate frequency‐dependent selection into the prey evolution, assuming that an individual's susceptibility to predation depends on both the trait it possesses and the traits of others in the population. When frequency‐dependent selection is symmetric, we show that the trait equation is unable to track changes in the fitness landscape, that is, the fitness landscape may change while the trait continues to evolve to the same trait value. As a result, evolution may produce unfit prey populations. Meanwhile, we show that asymmetric frequency‐dependence may have a destabilizing effect on the system, resulting in a closed invariant curve via a Neimark–Sacker bifurcation.

• Evolution of toxicant resistance may change system dynamics in a number of ways, such as enabling species persistence at higher toxicant levels. However, here we show that frequency-dependent evolution has the potential to destabilize system dynamics, resulting in stable closed invariant curves. • Destabilization of system dynamics may have important consequences for species persistence, as periods of low abundance may increase the possibility of species extinction.

K E Y W O R D S
discrete-time predator-prey models, Darwinian evolutionary game theory, frequency-dependent selection, Neimark-Sacker bifurcation 1 | INTRODUCTION Frequency-dependent selection occurs when an individual's fitness is dependent on the relative frequencies of the different phenotypes present in the population. Positive frequencydependent selection, which refers to the most common phenotypes having the greatest fitness, has been shown to lead to the evolution of warning coloration and mimicry (Takahashi & Kawata, 2013). Meanwhile, negative frequency-dependent selection, which occurs when rare phenotypes have a fitness advantage, has been credited with maintaining diversity. For instance, negative frequency-dependence was shown to extend the period of coexistence between two competing clades of Escherichia coli (Maddamsetti et al., 2015) and to maintain genetic polymorphism in the female-dimorphic damselfly Ischnura senegalensis (Takahashi & Kawata, 2013). Mathematically, it has been demonstrated that a stable polymorphic equilibrium is only possible if rare genotypes have a fitness advantage (Rueffler et al., 2006). In addition, it has been shown that asymmetric frequency-dependence, in which individuals with a trait larger than the mean trait experience reduced competitive pressure, may contribute to the evolution of semelparity (Cushing, 2019a). A number of mechanisms may lead to frequency-dependent evolution, including mating behavior, predator foraging strategies, mimicry, co-evolutionary interactions between populations, and niche heterogeneity (Udovic, 1980). Previously, the authors studied the dynamics of a discrete-time predator-prey model and an evolutionary counterpart (Ackleh et al., 2019;Ackleh, Hossain, et al., 2020). The interior equilibrium of this system was shown to be stable when it exists and was even shown to be globally asymptotically stable given certain conditions on the nullclines of the model (Ackleh, Salceanu, et al., 2020). The evolutionary model was used to examine how evolution in the prey species in response to an environmental disturbance or toxicant may impact the predator-prey dynamics. For this model, it was shown that, when the speed of evolution is sufficiently slow, the evolutionary dynamics can be described by a continuous perturbation of the nonevolutionary system. Meanwhile, fast evolution may destabilize the system dynamics resulting in cycles and even chaos.
A simplifying assumption of the evolutionary model described above is that an individual's fitness was assumed to only be impacted by the trait it possesses and not by the traits of others in the population. That is, this model did not consider how frequency-dependent evolution in the prey species may impact the dynamics of the predator-prey system. In this paper, we incorporate frequency-dependent evolution into the predator-prey system by assuming that a prey individual's susceptibility to predation contains a frequency-dependent component. For instance, this frequency-dependent predation may arise if the predator species exhibits frequency-dependent predation in which more common phenotypes are preferred. This behavior has been documented in a number of species, including the parasitoid wasp Anisopteromalus calandrae, which have demonstrated an increased preference for the more common host species as determined by olfactory cues (Ishii & Shimada, 2012), and blue jays (Cyanocitta cristata) which exhibit a preference for the more common type of cryptic Catocala moth (Bond & Kamil, 1998).
Here we consider two types of negative frequency-dependent selection. We show that, if this frequency-dependent component is maximized when an individual's trait is equal to the mean population trait (symmetric frequency-dependence), then the system dynamics are equivalent to the frequency-independent case studied in Ackleh et al. (2019) and Ackleh, Hossain, et al. (2020). However, this frequency-dependent component impacts the fitness landscape and may lead to the evolution of a less fit prey population. Alternatively, when the frequencydependent component describes asymmetric frequency-dependence such that individuals with a trait value larger than the mean trait experience enhanced predation risk, frequencydependent evolution of a prey species may destabilize the dynamics of a predator-prey system. This destabilization occurs via a different mechanism than when the speed of evolution is fast. In the latter, case destabilization follows a period-doubling-route-to-chaos in which variations in the population dynamics are small relative to variations in the mean trait. Meanwhile, with asymmetric frequency-dependence, destabilization may occur as a result of a Neimark-Sacker bifurcation. For a specific application, we show numerically that this bifurcation is supercritical, resulting in a stable closed invariant curve with significant variation in both the trait and population dynamics. This paper is organized as follows. In Section 2, we introduce the evolutionary predator-prey model with frequency-dependent selection. We then show in Section 3 that, on the boundary of the positive cone, the dynamics of this model are equivalent to the frequencyindependent case. We also show that, as with the previous model (Ackleh et al., 2019), the stability of the interior equilibrium may be obtained by taking the speed of evolution sufficiently slow. In Section 3.3 we discuss the conditions for a Neimark-Sacker bifurcation in the predator-prey model, and, in the appendix, show that the eigenvalue assignment conditions for a Neimark-Sacker bifurcation may be satisfied when frequency-dependence is asymmetric. Our theoretical results are corroborated in Section 4, where we use numerical simulation to show that a Neimark-Sacker bifurcation can, in fact, occur and results in the existence of a stable closed invariant curve. We also demonstrate how symmetric frequency-dependence may result in less fit prey populations since the trait evolution may be unable to respond to changes in the fitness landscape. Finally, in Section 5, we provide some concluding remarks.
where n and p denote the densities of the prey and predator, respectively. The unit of time is taken to be proportional to the generation time of the prey species, where the proportionality constant may be chosen based on the biology of a specific application. We assume that the prey population is able to evolve to resist a toxicant, where u represents the mean phenotypic trait possessed by the prey population and v denotes the trait of a focal individual. The traits of individuals in the prey population are assumed to be symmetrically distributed around the mean trait with constant variance ν > 0. This constant can be interpreted as defining the speed of evolution (Vincent & Brown, 2005). We define the prey growth rate as where r is the inherent prey growth rate, ϵ describes the reduction in the growth rate caused by the toxicant, ϕ gives the normalized (ϕ (0) = 1) density-dependent self-regulation of the prey species, and f is the probability that a prey individual is consumed by a predator individual. The predator equation assumes that predator reproduction is dependent on prey consumption, where b describes the conversion of consumed prey into new predator individuals. In addition, predator survival, s 0 < < 1, is assumed to be density-independent. Meanwhile, changes in the mean trait u are proportional to the fitness gradient of a focal individual, where fitness is taken to be R ln (Vincent & Brown, 2005). We assume that the nonlinearities ϕ and b are smooth functions in the following set X : This set includes Beverton-Holt-type nonlinearities. For further discussion on these model assumptions see Ackleh et al. (2019) and Ackleh, Hossain, et al. (2020).
To describe the evolutionary dynamics of model (1) we take the following trait dependencies: where r w w , , > 0 r 0 ϵ and 0 < ϵ < 1 0 . The first two nonlinearities assume that a trade-off occurs between the inherent growth rate and increased toxicant resistance such that both r and ϵ are maximized at trait value u = 0. Hence a trait value of u = 0 means that, on average, the population does not evolve resistance to the toxicant.
Unlike the model studied in Ackleh et al. (2019) and Ackleh, Hossain, et al. (2020), here we assume that the susceptibility of prey to predation, as defined by c v u ( , ), is frequency dependent. From a biological perspective, we may expect a prey individual's susceptibility to predation to depend both on the trait it possesses, as well as how that trait compares to the traits of others. For example, evolution of resistance to cadmium in Daphnia magna reduces body size, which in turn decreases swim speed (Ward & Robinson, 2005). This means an individual that has evolved greater toxicant resistance than other prey individuals may have an enhanced risk of predation since it is slower than the others (i.e., an individual's predation risk depends on its trait relative to others). However, certain predators of Daphnia, such as the phantom midge Chaoborus, are limited by gape size and cannot eat prey that are too large (Riessen & Trevett-Smith, 2009). Thus, predation risk also depends on the individual trait (in this case size). Motivated by this, here we define the susceptibility of prey to predation as

| Existence and stability of boundary equilibria
We first consider the evolution of the mean trait when it is independent of both predator and prey densities. Since the trait equation for u is independent of n, this may be obtained by setting p = 0. Notice that, if p = 0, then the trait value is determined by the trade-off between the growth rate and toxicant effect, that is, evolution is independent of c and we have the same dynamic outcomes for both cases. When p = 0, the trait equation is given by In what follows, we restrict our focus to nonnegative trait values. Notice that, from Equation (6) Manipulation of this equation shows that this inequality is satisfied when When inequality (7) holds, Equation (6) has two equilibria, u = 0 and the positive equilibrium In Lemma 3.1, we show that the positive equilibrium is globally asymptotically stable when it exists.
, then Equation (6) has exactly one equilibrium u = 0 which is globally Then, in addition to the equilibrium u = 0, Equation (6) has a unique positive equilibrium ũ. The equilibrium u = 0 is unstable and there exists a ν* > 0 such that the equilibrium ũ is globally asymptotically stable for ν ν 0 < < *. Proof.
where g is the map given by the right-hand side of (6). Then | | g′(0) < 1, that is, u = 0 is locally asymptotically stable, when Therefore, g u u ( ) < for any u > 0. From Theorem 2.5 in Allen (2007), it follows that u = 0 is globally asymptotically stable.
. Then the positive equilibrium ũ exists and | | g u ′(˜) < 1, that is, ũ is locally asymptotically stable, when >˜, and g does not have a maximum on u (0,˜). By Theorem 2.6 in Allen (2007), ũ is globally asymptotically stable if and only if Equation (6) has no 2-cycles. Since ≠ g u νh u 1 + ′( ) = 2 + ′( ) 0 for ν ν < *, by Theorem 2.7 in Allen (2007), Equation (6) has no 2-cycles. Thus, ũ is globally asymptotically stable. □ Next we consider the existence and stability of the boundary equilibria. We again note that the symmetric and asymmetric frequency-dependent cases are the same except when ≠ p 0.
and is locally asymptotically stable if w w w and s b n nc ϵ < + + (¯)¯< 1.
, and is locally asymptotically stable if ( ) ν w w w w and r u u < 2 ln ( + ) and r u u (˜)(1 − ϵ(˜)) > 1, and is locally asymptotically stable if ( ) ν w w w w and s b n nc u < 2 ln ( + ) Remark 4. The proof of Theorem 3.2 follows from straightforward linearization and eigenvalue calculations. Since this result is analogous to Theorem 3.3 in Ackleh et al.

| Interior dynamics
Given the high-order nature of the equations governing the interior dynamics, it is difficult to study the interior dynamics through direct computation. Instead, we study the interior dynamics via a combination of persistence theory, perturbation analysis, and numerical simulation.
To study the interior dynamics, we use the notation or >0 for c given by (3a) or (3b), respectively. First, we show that there is a compact positively invariant set that attracts solutions to model (1), that is, there exist positive constants N P , , and U such that all solutions starting with nonnegative initial conditions eventually enter the compact set To this end, we write the trait equation as u t k u ν u t From this we can argue as in In Lemma 3.3 and Theorem 3.4, we establish conditions for the persistence of system (1), that is, we show that solutions are bounded away from the boundary of the positive cone. These results are again analogous to those developed in Ackleh et al. (2019).
for any positive initial condition.
, there exist positive ε and σ such that for any u (0) > 0. By Proposition 3.2 in Magal and Zhao (2005) Then model (1) is persistent, that is, there exists a positive ε such that for any positive initial condition.
Once the persistence of the trait u is established, the persistence of the prey and predator follows the same arguments as used in the proof of Theorem 3.12 in Ackleh et al. (2019). Therefore, we omit the proof of Theorem 3.4 here.
A consequence of the existence of a compact positively invariant attracting set and Theorem 3.4 is that, if the inequalities in (11) hold then system (1) is permanent and has an interior equilibrium (see Hutson, 1990). This equilibrium is stable for sufficiently small ν, as is shown in Theorem 3.5. Though the proof of this theorem is similar to the proof of Lemma 3.7 in Ackleh et al. (2019), we give the details in full in Appendix A as the notations introduced in this proof will be used to study the existence of a Neimark-Sacker bifurcation.

| Neimark-Sacker bifurcation
A Neimark-Sacker bifurcation occurs when a pair of complex eigenvalues leave the unit circle. However, given that the eigenvalues for the Jacobian matrix associated with the interior equilibrium are the solutions to a cubic polynomial, it is not convenient to work with the exact expressions for these eigenvalues. Instead, to study the existence of a Neimark-Sacker bifurcation we apply a result given in Wen (2005). For convenience, this result is provided in Theorem B.1 in Appendix B.
In Appendix B, we show in Theorem B.2 that the four eigenvalue assignment conditions listed in Theorem B.1 may be satisfied by fixing ν sufficiently small and increasing the value of w c . Thus, Theorem B.2 suggests that it may be possible to obtain a Neimark-Sacker bifurcation by increasing the value of w c from zero, that is, changing evolution in model (1) from frequency-independent selection to asymmetric frequency-dependent selection. Motivated by this, we take the bifurcation parameter to be w c and the control parameter, as defined in Theorem B.1, to be ν. Using the notation that was introduced in the proof Theorem 3.5, the conditions given in Theorem B.1 for model (1)  (1) > 0, (−1) < 0, In Section 4, we use numerical simulation to further study the existence of a Neimark-Sacker bifurcation.

| APPLICATIONS
In this section, we use numerical simulations to explore the implications of symmetric and asymmetric frequency-dependence. For all numerical simulations, we define ϕ and b according to the Beverton-Holt nonlinearities

| Symmetric frequency-dependence
As we noted in Section 2, when model (1) has symmetric frequency-dependence, as defined by Equation (3a), both the population and trait dynamics are independent of the parameter w c (see Equations 4 and 5a). However, the adaptive landscape R n p v u ( , , , ) e e e for a given equilibrium n p u ( , , ) e e e does depend on w c . As a consequence, as the parameter w c is varied, the equilibria and their stability do not change but the critical points of the adaptive landscape may change. Since evolution is unable to track these changes in the adaptive landscape, it is possible for symmetric frequency-dependence to produce unfit prey populations, that is, populations whose mean trait is at a fitness minimum of the adaptive landscape. This is illustrated in Example 4.1.
Example 4.1 (Symmetric frequency-dependence may produce unfit populations). Figure 1a and 1b gives the time-series dynamics for the predator and prey densities and mean trait, respectively, for the parameter values:  Figure 1c gives the adaptive landscape for these parameter values along with two values of w c : w = 0.1 c and 0.9. We can see that, as w c is increased, the fitness maximum becomes a fitness minimum, resulting in a population whose fitness is not maximized.
Recall that changes in the mean trait are determined by the fitness gradient, Mathematically, we observe that the trait equation is independent of the parameter w c because the assumption that the frequency-dependent component c v u ( , ) f is maximized at the mean trait results in its partial derivative vanishing when evaluated at the mean trait, However, the trait equation was obtained by approximating the fitness of an individual with trait u j by a first-order Taylor expansion around the mean trait u. Therefore, if we instead approximate this term using a second-order expansion, then the trait equation is no longer independent of w c . Specifically, the trait equation becomes see Appendix C for details on the derivation of this equation. In particular, we note that the assumption the traits are symmetrically distributed around the mean trait u means that the additional term in this equation is of the same order with respect to the speed of evolution ν.

n t p t v u t R n t p t v u t u t R n t p t v u t R n t p t v u t
With the addition of the third term in trait equation (13), the traditional interpretation of an equilibrium trait value occurring at a critical point of the fitness landscape no longer holds. We also observe a new dynamic, often termed evolutionary suicide (Gyllenberg & Parvinen, 2001;Nonaka et al., 2013), in which evolution results in population extinction. Such a dynamic cannot occur when symmetric frequency-dependence is given by the traditional trait equation. This is illustrated in Example 4.2.
Example 4.2 (Symmetric frequency-dependence with a modified trait equation may lead to evolutionary suicide). Here we consider the new trait equation (13) with the same parameter values used in Example 4.1. Figure 2a gives the fitness landscapes and trait equilibrium values for four values of w c . We now observe that, when w c is small, a trait equilibrium value corresponds to a point that is close to but below the maximum fitness (compare w = 0.1 c to Figure 1c). As w c increases, the mean trait approaches the fitness maximum while the fitness landscape starts to flatten out. Eventually, for w = 0.65 c , we observe that the mean trait overshoots the fitness maximum. For even larger values of w c , F I G U R E 2 (a) The fitness landscapes for fours values of w c when evolution is defined by Equation (13) which correspond to fitness minimums when using trait equation (1c), the trait component does not equilibrate. Instead, the mean trait value continues to increase, eventually resulting in both populations crashing, as shown in Figure 2b and 2c.

| Asymmetric frequency-dependence
In Theorem B.2, we showed that the eigenvalue assignment conditions for a Neimark-Sacker bifurcation may hold for the asymmetric frequency-dependence case, using w c as the bifurcation parameter. In Example 4.3, we provide an example in which increasing the parameter w c leads to a supercritical Neimark-Sacker bifurcation resulting in a stable closed invariant curve. In Example 4.4, we contrast these oscillations with the 2-cycles that may be obtained when the speed of evolution ν is fast.  Finally, we follow Kuznetsov (2013) and Xin et al. (2010) to determine that the direction of bifurcation is supercritical since d, as defined in Appendix B, is given as d = −0.14694 < 0. Thus a supercritical Neimark-Sacker bifurcation occurs at w c 0 resulting in a stable closed invariant curve. This is illustrated in Figure 3 which gives the time series solution for model (1) with the parameter values listed above for two values of w c . For w = 0.14 c , it is shown in Figure 3a and 3b that the model exhibits a stable interior equilibrium. However, if we increase w c to w = 0.5 c , this interior equilibrium destabilizes resulting in quasiperiodic oscillations, as shown in Figure 3c and 3d. To better illustrate these oscillations, in Figure 4a we give the invariant curve in phase space obtained from plotting the last 500 values of a 5000 iteration simulation. Finally, in Figure 4b we plot the absolute value of the eigenvalues of the Jacobian matrix evaluated at the interior equilibrium as a function Observe that a pair of complex eigenvalues leaves the unit circle at w c 0 , after which no eigenvalues enter or leave the unit circle. This suggests that the interior equilibrium does not experience any additional bifurcations.
Example 4.3 shows that asymmetric frequency-dependence produces new mathematical dynamics not observed in the symmetric case. However, from an ecological perspective, the quasiperiodic oscillations observed here could be thought of as another form of evolutionary suicide. Specifically, these oscillations result in periods of time over which the prey population is at low density. Should some type of stressor be applied to the population during these times, this could result in the prey population dying out. For comparison purposes, in Figure 5 we consider how model dynamics are changed in the asymmetric frequency-dependence case when the trait evolution is described by the modified trait equation (13). Figure 5 gives the time-series dynamics for the densities and mean trait when trait evolution is described by Equation (13) with the same parameter values used in Figure 3c and 3d. Simulation shows that stable quasiperiodic oscillations also occur in this case. As is evident from Figure 5, the modified trait equation increases the possibility of evolutionary suicide as it results in extended periods of time both species spend at low densities. Previously, we showed that the nonevolutionary version of model (1) (which can be obtained by fixing ≡ u t u ( )¯and setting ν = 0) does not possess 2-cycles (Ackleh, Hossain, et al., 2020). By continuity, this holds for the evolutionary model provided that ≈ ν 0. However, we showed numerically that increasing ν may lead to a bifurcation resulting in stable 2-cycles. In Example 4.4, we demonstrate how the Neimark-Sacker bifurcation that occurs with asymmetric frequency-dependence presents, from an ecological perspective, a fundamentally different type of behavior when compared with the 2-cycles that result from fast evolution. While both scenarios may be interpreted as "boom-bust" dynamics, the Neimark-Sacker bifurcation results in longer periods of time at lower densities, leading to an enhanced risk of prey extinction due to stochastic effects.
Example 4.4 (Stable 2-cycles resulting from fast evolution). In Figure 6, we provide two examples of stable 2-cycles that result from increasing the speed of evolution ν. In Figure 6a and 6b, we first consider the case when all parameter values except ν are the same  Figure 3a and 3b. As we showed in Example 4.3, when ν is sufficiently small, these parameters result in a stable interior equilibrium. However, Figure 6 shows that increasing ν results in a stable 2-cycle. Notice that the variation in the population densities is much smaller when compared with the oscillations arising from a Neimark-Sacker bifurcation. In addition, the densities do not reach low levels that may indicate an extinction risk due to stochastic effects. In Figure 6c and 6d we present another example in which greater variation is observed in both the trait and population densities, but this variation would still be unlikely to result in an extinction event. Even larger amplitudes in the prey 2-cycle are possible if parameter values are chosen so that the predator density is decreased (not shown). However, these cycles still do not bring the prey population down to low densities.

| CONCLUSION
In Ackleh et al. (2019) we proved that, for the nonevolutionary version of model (1) (obtained by fixing u as a parameter), the interior equilibrium is stable when it exists. However, we showed in Ackleh, Hossain, et al. (2020) that (frequency-independent) evolution may cause the interior equilibrium to destabilize when the speed of evolution, as described by ν, is fast. This destabilization results in a period-doubling-route-to-chaos in which the mean trait exhibits significant variation but this results in only minor oscillations in the predator and prey populations. Here we show that frequency-dependent evolution may also destabilize the population dynamics. Unlike the previous case, this destabilization may occur for small ν and is the result of two complex eigenvalues, rather than a single real eigenvalue, leaving the unit circle. Thus, frequency-dependent evolution introduces new dynamics not observed in either the nonevolutionary or the frequency-independent evolutionary model. Moreover, as demonstrated in Example 4.3, we see that frequency-dependent evolution may introduce significant oscillations in both the trait and population densities. These oscillations may be detrimental to both the populations since they result in extended periods of time at which the populations are at low densities. We also observe that, when frequency-dependent selection is symmetric, the trait equation is unable to track changes in the fitness landscape, which may result in the prey population evolving to a fitness minimum, as is shown in Example 4.1. This disconnect between trait evolution and fitness results from approximating the fitness of an individual with a given phenotype by a first-order expansion around the mean trait. Motivated by this, we developed a modified trait equation using a second-order approximation so that trait evolution is no longer independent of the frequency-dependent parameter w c . In Example 4.2 we observe that the second-order approximation produces new behavior that is not captured with the first-order approximation. Specifically, Figure 2 shows that adaptation may lead to a population crash.
The idea that adaptation of the individual may sometimes lead to selection of traits having catastrophic results on the population was suggested by Haldane in 1932(Haldane, 1932. Such a phenomenon, known as evolutionary suicide (Gyllenberg & Parvinen, 2001;Nonaka et al., 2013), has been theoretically demonstrated to occur using mathematical models (Gyllenberg & Parvinen, 2001;Nonaka et al., 2013) but several recent empirical studies have also suggested that this phenomenon may occur in nature (Rankin & Lopez-Sepulcre, 2005). A similar phenomenon to what is observed in Example 4.2, that is, population crashing as a result of an increasing trait, was termed evolutionary suicide as a result of 'runaway selection' in Nonaka et al. (2013). A notable distinction between the two types of frequency-dependent selection studied here is that, while the oscillations observed in the asymmetric frequency-dependent case may also be interpreted as a form of evolutionary suicide, this phenomenon cannot occur in the symmetricfrequency-dependent case when evolution is described by a first-order expansion. Notice that for ν δ < 1 the only quantity that contributes positively to D (−1) is J J J −˜˜1 3 21 32 . However, since this term is linear in ν, there exists a δ > 0 2 such that D (−1) < 0 for ν δ δ < min{ , } 1 2 . To consider the last Jury condition, we first introduce the simplifying notation:   (1) with asymmetric frequency-dependence. We also review how to determine the direction of bifurcation.
Theorem B.1 (Wen, 2005). Consider the n-dimensional system , where ∈ μ  is a bifurcation parameter. Assume that this system has an equilibrium x* and the characteristic polynomial of the Jacobian matrix evaluated at x* is given by where a a μ K i n = ( , ), = 1, …, i i for some control parameter ∈ K . Define the sequence of and assume that the following conditions hold: (i) Eigenvalue assignment: Then a Neimark-Sacker bifurcation occurs at μ μ = 0 .
In Example 4.3, we numerically calculate the quantity d which determines the direction of bifurcation for the observed Neimark-Sacker bifurcation. Here we follow Kuznetsov (2013) and Xin et al. (2010) to outline the steps required to calculate this quantity. We first make a change of variables so that the interior equilibrium n p u (ˆ,ˆ,ˆ) is transformed to the origin. This is can be accomplished using .