Chaotic dynamics in a simple predator-prey model with discrete delay
aa r X i v : . [ m a t h . D S ] J u l Chaotic dynamics in a simple predator-prey model with discretedelay
Guihong Fan a and Gail S. K. Wolkowicz b a Department of Mathematics, Columbus State University, Columbus, Georgia 31907; b Department of Mathematics and Statistics, McMaster University, Hamilton, Ontario,Canada L8S 4K1
ARTICLE HISTORY
Compiled August 3, 2020
ABSTRACT
A discrete delay is included to model the time between the capture of the prey andits conversion to viable biomass in the simplest classical Gause type predator-preymodel that has equilibrium dynamics without delay. As the delay increases fromzero, the coexistence equilibrium undergoes a supercritical Hopf bifurcation, twosaddle-node bifurcations of limit cycles, and a cascade of period doublings, eventu-ally leading to chaos. The resulting periodic orbits and the strange attractor resembletheir counterparts for the Mackey-Glass equation. Due to the global stability of thesystem without delay, this complicated dynamics can be solely attributed to theintroduction of the delay. Since many models include predator-prey like interactionsas submodels, this study emphasizes the importance of understanding the impli-cations of overlooking delay in such models on the reliability of the model-basedpredictions, especially since temperature is known to have an effect on the length ofcertain delays.
KEYWORDS
Predator-prey model; stage-structured model with maturation delay; Hopf andsaddle-node bifurcation of limit cycles; Period doubling route to chaos; bi-stability;Mackey-Glass attractor; uniform persistence
AMS (MOS) subject classification:
1. Introduction
A Gause type predator-prey model with response function f ( x ) is given by ˙ x ( t ) = rx ( t ) (cid:18) − x ( t ) K (cid:19) − y ( t ) f ( x ( t )) , ˙ y ( t ) = − sy ( t ) + Y y ( t ) f ( x ( t )) , (1)where x ( t ) denotes the density of the prey population and y ( t ) the density of predators.Parameters r, K, s , and Y are positive constants denoting the intrinsic growth rate CONTACT G. S. K. Wolkowicz. Email: [email protected] nd the carrying capacity of the prey, the death rate of the predator in the absenceof prey, and the growth yield constant for the conversion of prey to viable predatordensity, respectively.If f ( x ) is of Holling type I form in model (1) (i.e. f ( x ) = mx where m is a positiveconstant denoting the maximal growth rate of the predator), it is well-known (seee.g. [2, 9]) that either the predator population approaches extinction and the preypopulation approaches its carrying capacity, or the predator population and the preypopulation coexist and their density approaches a positive equilibrium. Hence, for allchoices of the parameters, all solutions of this system approach a globally asymptot-ically stable equilibrium, and so any nontrivial oscillatory behaviour that arises dueto the introduction of delay in the model can be attributed solely to the delay. Forthis reason, we choose Holling type I response functions instead of the more realisticHolling type II form, since the Holling type II form results in a model that gives rise tonontrivial period solutions without delay (see Rosenzweig [25]). One would also expectthat any exotic dynamics that the model with Holling type I form admits due to theintroduction of delay would be shared by the model with Holling type II form. Li etal. [19] studied this model with the Holling type II response function of Monod formand showed that stability switches caused by varying the time delay are accompaniedby bounded global Hopf branches, and they proved that when multiple Hopf branchesexist, they are nested and the overlap produces coexistence of two or possibly morestable limit cycles. However, they did not go on to discover the even richer dynamicsthat our analysis suggests exists in that case.Incorporating a time delay in (1) to model the time between the capture of the preyby the predator and its conversion to viable predator biomass, in the case of Hollingtype I functional response, f ( x ) = mx, m >
0, we obtain the following system: ˙ x ( t ) = rx ( t ) (cid:18) − x ( t ) K (cid:19) − my ( t ) x ( t ) , ˙ y ( t ) = − sy ( t ) + Y e − sτ my ( t − τ ) x ( t − τ ) . (2)The term e − sτ y ( t − τ ) represents those predators that survive the τ ≥ t − τ in the past. Thus, we haveincorporated the delay in the growth term of the predator equation in a manner thatis consistent with its decline rate given by the model, as described in Arino, Wang,and Wolkowicz [1].We define R + ≡ { x ∈ R : x > } , int R + ≡ { x ∈ R : x > } . We denote by C ([ − τ, , R + ), the Banach space of continuous functions from the interval [ − τ,
0] into R + , equipped with the uniform norm. We assume initial data for model (2) is takenfrom X = C ([ − τ, , R + ) × C ([ − τ, , R + ) . (3)Model (2) also has other interpretations. Gourley and Kuang [12] studied a stage-structured predator-prey model in which they included an equation for the juvenilepredators and assumed a constant maturation time delay, i.e., they assumed thatthe juvenile predators take a fixed time to mature. Using the approach developed inBeretta and Kuang [3], the authors considered the possibility of stability switches, andconcluded that there is a range of the parameter modeling the time delay for whichthere are periodic solutions. If the juveniles in their model suffer the same mortality2ate as adult predators, their model decouples and yields model (2). Forde [8] alsoconsidered this model and conjectured that there are periodic orbits whenever theinterior equilibrium exists and is unstable. He also noted that if the interior equilibriumexists and is asymptotically stable without delay, then for small delays it remainsglobally asymptotically stable. We show that whenever the interior equilibrium existsand is unstable, the system is uniformly persistent. Gourley and Kuang [12] had alreadyshowed that a Hopf bifurcation eventually occurs if the delay is increased, destabilizingthis equilibrium and giving birth to a nontrivial periodic solution. Forde [8] left asan open question whether more than one periodic orbit is possible and provided anumerical example suggesting chaos is possible but does not consider the route tochaos. We give numerical evidence that there is a range of parameters for which twostable periodic orbits and an unstable periodic orbit all exist and we show a period-doubling route to chaos followed by a period-halfing route back to stability of theinterior equilibrium.Cooke, Elderkin, and Huang [4] considered a model similar to the one in Gourley andKuang [12], and obtained results concerning Hopf bifurcation of a scaled version. Thescaling they used eliminated the parameter modelling the time delay, the parameterthat we focus on and use as a bifurcation parameter. This simplified their analysis, sincethen, unlike in our case, the components of the coexistence equilibrium are independentof the time delay.In this manuscript we show that the introduction of time delay cannot only desta-bilize the globally asymptotically stable coexistence equilibrium of model (2), it canalso be responsible for exotic dynamics for intermediate values of the delay as wellas the eventual disappearance of the coexistence equilibrium with the extinction ofthe predator for large enough delays. Although chaotic dynamics has been observedin other models of predator-prey interactions, the other models either require at leastthree trophic levels, or the response functions are not as simple and so the models ad-mit oscillatory behavior even in the absence of delay, or the other models incorporatethe delay in such a way that the predators still contribute to population growth evenif the time required to process the prey is longer than the life-span of the predator(i.e., the factor e − sτ is missing in the ˙ y equation), or the delay is used to model differ-ent mechanisms (see e.g., [11, 15, 24, 28]). The observation that the resulting strangeattractor resembles the strange attractor for the Mackey-Glass equation [21] is alsonew.This paper is organized as follows. In section 2, we scale the model and show thatit is well-posed. In section 3, we consider the existence and stability of equilibria. Ifparameters are set so that it is possible for the predator to survive when there is nodelay, it is well-known that the equilibrium at which both the prey and the predatorsurvive is globally asymptotically stable with respect to positive initial conditions(i.e. solutions approach this equilibrium for any choice of positive initial data). In thecase of delay, the components of this coexistence equilibrium depend on the delay.We prove that for positive delay, when this equilibrium exists, both the predator andthe prey populations persist uniformly. However, a sufficiently long delay results inthe disappearance of this equilibrium, resulting in the extinction of the predator andconvergence to a globally asymptotically stable equilibrium with the prey at carryingcapacity. We give criteria which when satisfied imply that there are at least two Hopfbifurcations that occur before the extinction of the predator, resulting in sustainedoscillatory behaviour for intermediate values of the delay. Finally, in section 4, bymeans of time series, time delay embeddings, and orbit (bifurcation) diagrams weshow that there are saddle-node bifurcations of limit cycles resulting in bistability3s well as sequences of period doubling bifurcations leading to chaos, with a strangeattractor resembling the strange attractor for the Mackey-Glass equation [21]. Weconclude with a brief discussion.
2. Scaling and basic properties of solutions
In order to simplify the analysis, we introduce the following change of variables:˘ t = rt, ˘ x (˘ t ) = x ( t ) /K, ˘ y (˘ t ) = my ( t ) /r, ˘ τ = rτ, ˘ s = sr , ˘ Y = Y Km/r. (4)We drop the˘’s for convenience and study the equivalent scaled version of model (2): ˙ x ( t ) = x ( t )(1 − x ( t )) − y ( t ) x ( t ) , ˙ y ( t ) = − sy ( t ) + Y e − sτ y ( t − τ ) x ( t − τ ) , ( x ( t ) , y ( t )) = ( φ ( t ) , ψ ( t )) ∈ X, for t ∈ [ − τ, , (5)where X was defined in (3).First we address well-posedness of system (5). For positive delay τ , the existenceand uniqueness of solutions of system (5) was shown in Gourley and Kuang [12]. Thefollowing proposition, proved in A.1, indicates that for positive delay the solutionsremain nonnegative and provides an upper bound for each component. Proposition 2.1.
Consider model (5) with initial data in X .(1) The solutions exist, are unique, and remain nonnegative for all t > .(2) lim sup t →∞ x ( t ) and lim sup t →∞ y ( t ) s Y e − sτ ( s + 1) .(3) Consider model (5) with initial data in X where X = { ( φ ( t ) , ψ ( t )) ∈ X : φ (0) > ∃ θ ∈ [ − τ, s.t. φ ( θ ) ψ ( θ ) > . } (6) Then, x ( t ) > for all t > and there exists T > such that y ( t ) > for all t > T .
3. Existence and stability of equilibria and uniform persistence
Model (5) can have up to three distinct equilibria: E = (0 , , E = (1 , , E + = ( x + ( τ ) , y + ( τ )) = (cid:16) sY e sτ , − sY e sτ (cid:17) . (7)The components of E + are nonnegative and E + is distinct from E , if, and only if,0 τ < τ c , where τ c = 1 s ln (cid:18) Ys (cid:19) . (8)4hus, when τ c >
0, i.e., when
Y > s , the components of E + are both positive, and E + is referred to as the coexistence equilibrium.When there is no delay, i.e. τ = 0 in (5), E is always a saddle attracting solutionswith x (0) = 0. If Y < s , one of the components of E + is negative and so it is notrelevant, and E is globally asymptotically stable with respect to initial conditionssatisfying x (0) > y (0) >
0. When Y = s , E and E + coalesce and are globallyattracting provided x (0) >
0. If
Y > s , then E is a saddle attracting solutions with x (0) > y (0) = 0 and E + sits in int R and is global asymptotically stable withrespect to initial conditions in int R .When τ >
0, to determine the local stability of each equilibrium solution, we usethe linearization technique for differential equations with discrete delays (see Haleand Lunel [13]). After linearizing (5) about any one of these equilibria, ( x ⋆ , y ⋆ ), thecharacteristic equation, P ( λ ) | ( x ⋆ ,y ⋆ ) = 0, is given by,( λ + s )( λ + y ⋆ − (1 − x ⋆ )) + Y e − ( s + λ ) τ x ⋆ (1 − x ⋆ ) − λY e − ( s + λ ) τ x ⋆ = 0 . (9)We summarize the results on local and global stability of the equilibrium pointsand uniform persistence of the populations in the following theorem. The proof canbe found in A.2. Theorem 3.1.
Consider (5).(1) Equilibrium E is always unstable.(2) Equilibrium E is(a) unstable if τ < τ c , and(b) globally asymptotically stable (with respect to X ) if τ > τ c , (i.e. if se sτ Y > ).(3) Both components of E + are positive (i.e., E + exists), if, and only if, τ < τ c ,(i.e. se sτ Y < ).(a) When E + exists and τ = 0 , E + is globally asymptotically stable with respectto int R .(b) When E + exists, and τ > , model (5) is uniformly persistent with respectto initial data in X , i.e., there exists ǫ > independent of ( φ ( t ) , ψ ( t )) ∈ X such that lim inf t →∞ x ( t ) > ǫ and lim inf t →∞ y ( t ) > ǫ . Thus, for any fixed time delay τ , if sY e − sτ >
1, only the prey population survivesand it converges to a steady state. On the other hand, if the inequality is reversed,for appropriate initial data both the prey and the predator populations are uniformlypersistent, i.e. survive indefinitely. However, we have not yet addressed what form thedynamics takes in the latter case. E + When E + exists, by Theorem 3.1 both populations survive indefinitely. To address thepossible forms the dynamics can take, we begin by investigating the local stability of E + when it exists, i.e., when 0 τ < τ c , and hence both components are positive.Evaluating the characteristic equation (9) at E + gives P ( λ ) | E + = λ + λs (cid:18) e sτ Y (cid:19) + s Y e sτ + e − λτ s (cid:18) − λ + (cid:18) − se sτ Y (cid:19)(cid:19) = 0 . P ( λ ) | E + = 0 is of the form P ( λ ) | E + = λ + p ( τ ) λ + ( qλ + c ( τ )) e − λτ + α ( τ ) = 0 , (10)where p ( τ ) = s (cid:18) e sτ Y (cid:19) , q = − s, c ( τ ) = s (cid:18) − se sτ Y (cid:19) , and α ( τ ) = s e sτ Y , (11)First assume that τ = 0. Then (10) reduces to λ + ( p (0) + q ) λ + ( α (0) + c (0)) = 0 . Since α (0) + c (0) = s (cid:0) − sY (cid:1) = sy + (0) > p (0) + q = se sτ Y >
0, by the Routh-Hurwitz criterion [10], all roots of (10) have negative real part. Therefore, E + is locallyasymptotically stable when τ = 0 and hence also for τ > E + as τ varies in the interval 0 < τ < τ c . Here, P (0) | E + = α ( τ ) + c ( τ ) = s y + ( τ ) > λ = 0 is not a root of (10). Therefore, theonly ways that E + can lose stability is: (i) when one of the characteristic roots equalszero. This only occurs when τ = τ c . This gives rise to a transcritical bifurcation where E + coalesces with E and then disappears as τ increases through τ c ; (ii) if characteristicroots bifurcate in from infinity; or (iii) if a pair of complex roots with negative realparts and non-zero imaginary parts cross the imaginary axis as τ increases from 0,potentially resulting in Hopf bifurcation. In A.3 we prove that (ii) is impossible toobtain the following lemma. Lemma 3.2. As τ increases from zero, the number of roots of (10) with positive realpart can change only if a root appears on or crosses the imaginary axis as τ varies. In order to determine when Hopf bifurcations occur, we first determine for whatvalues of τ pure imaginary roots of (10) exist so that (iii) can occur. We will also beinterested in secondary Hopf bifurcations.Suppose that λ = iω ( ω >
0) is a root of P ( λ ) | E + = 0, where i = √−
1. Then P ( iω ) | E + = − ω + ip ( τ ) ω + ( iqω + c ( τ )) e − iτω + α ( τ ) = 0 . Using Euler’s identity, e iθ = cos θ + i sin θ , and equating the real and imaginary parts,this is equivalent to c ( τ ) cos( τ ω ) + qω sin( τ ω ) = ω − α ( τ ) ,c ( τ ) sin( τ ω ) − qω cos( τ ω ) = p ( τ ) ω. Solving for cos( τ ω ) and sin( τ ω ) givessin( τ ω ) = c ( τ )( p ( τ ) ω ) + qω ( ω − α ( τ )) c ( τ ) + q ω , (12a)cos( τ ω ) = c ( τ )( ω − α ( τ )) + qω ( − p ( τ ) ω ) c ( τ ) + q ω . (12b)6quaring both sides of the equations in (12), adding, and rearranging gives ω + ( p ( τ ) − q − α ( τ )) ω + α ( τ ) − c ( τ ) = 0 . (13)Noting that (13) is a quadratic function of ω , we use the quadratic formula to obtain ω ± ( τ ) = 12 (cid:18) q − p ( τ ) + 2 α ( τ ) ± p ( q − p ( τ ) + 2 α ( τ )) − α ( τ ) − c ( τ )) (cid:19) . Substituting using (11), it follows that ω ± ( τ ) = 12 (cid:18) − (cid:18) se sτ Y (cid:19) ± s(cid:18) se sτ Y (cid:19) + s (cid:18) s e sτ Y − se sτ Y + 4 (cid:19)(cid:19) . (14)In order to determine for what values of τ there are positive real roots of (13), andhence candidates for pure imaginary roots, and possibly Hopf bifurcations, we define τ ∗ = 1 s ln (cid:18) Y s (cid:19) . (15)We will prove that for τ > τ ∗ , there are no postive real roots and for 0 τ < τ ∗ ω + ( τ ) = vuut (cid:18) − (cid:18) se sτ Y (cid:19) + s(cid:18) se sτ Y (cid:19) + s (cid:18) s e sτ Y − se sτ Y + 4 (cid:19)(cid:19) (16)is the only positive real root. Remark 1.
Note that τ ∗ >
0, if, and only if, sY < , and then x + ( τ ∗ ) = se sτ ∗ Y = .In the following theorem, proved in A.4, we give necessary conditions on τ for Hopfbifurcations to occur. Theorem 3.3.
Consider (5). Assume that both components of E + are positive.(1) ω + ( τ ∗ ) = 0 . If τ > τ ∗ , then (13) has no positive real root. Therefore, there canbe no pure imaginary roots of (10) in this case. In particular, if sY > , then τ ∗ , and so there can be no Hopf bifurcation of E + for any τ > .(2) Assume that sY < , and hence, τ ∗ > . If τ ∈ [0 , τ ∗ ) , then x + ( τ ) ∈ [ sY , ) , and(13) has exactly one positive real root, ω + ( τ ) , given by (16). If (10) has pureimaginary roots at τ , and hence τ is a candidate for Hopf bifurcation of E + ,then τ ∈ (0 , τ ∗ ) and ω + ( τ ) must satisfy (16). Remark 2. (1) Note that even though ω + (0) > sY < , τ = 0 is not acandidate for a Hopf bifurcation, since in this case all roots of (10) have negativereal parts. Also, τ ∗ is not a candidate, since ω + ( τ ∗ ) = 0.(2) Haque [14] finds an expression for ω + ( τ ∗ ) for a different scaling of the model.However, he does not go on as we do in what follows, to determine when theequations given in (12) are both simultaneously satisfied for ω + ( τ ∗ ), and to showthat the Hopf bifurcations are nested and how the number of Hopf bifurcationsincreases as the death rate of the predator decreases.7ubstituting the values of the coefficients given by (11), in the right-hand side of(12), and recalling that x + ( τ ) = se sτ /Y , we define h ( ω, τ ) = ωs (cid:18) s + x + ( τ ) − sx + ( τ ) − x ( τ ) − ω (1 − x + ( τ )) + ω (cid:19) , (17a) h ( ω, τ ) = ω (1 + s − x + ( τ )) − (1 − x + ( τ )) sx + ( τ ) s ((1 − x + ( τ )) + ω ) . (17b)Properties of the functions h and h are summarized in A.5. Remark 3.
By Theorem 3.3 and Lemma A.1, τ satisfies (12) for ω >
0, (and hence(10) has a pair of pure imaginary eigenvalues) if and only if τ ∈ (0 , τ ∗ ), andsin( τ ω + ( τ )) = h ( ω + ( τ ) , τ ) , (18a)cos( τ ω + ( τ )) = h ( ω + ( τ ) , τ ) . (18b)Define the function θ : [0 , τ ∗ ] → [0 , π ] (19) θ ( τ ) := arccos( h ( ω + ( τ ) , τ )) . By part 3 of Lemma A.1, stated and proved in A.5, θ ( τ ) is a well-defined, continuouslydifferentiable function. Replacing τ ω + ( τ ) by θ ( τ ) + 2 nπ in the left hand side of (18),we obtain sin( θ ( τ ) + 2 nπ ) = h ( ω + ( τ ) , τ ) , (20a)cos( θ ( τ ) + 2 nπ ) = h ( ω + ( τ ) , τ ) . (20b)Equation (20b) is satisfied directly by the definition of θ ( τ ). Equation (20a) is alsosatisfied, from parts 1 and 2 of Lemma A.1, since 0 θ ( τ ) π. By comparing (18) and (20), it follows that solutions of (18) occur at precisely thosepoints where the curves τ ω + ( τ ) and θ ( τ ) + 2 nπ intersect. Remark 4.
By Remark 3, (10) has a pair of pure imaginary roots at precisely thosepoints where the curves τ ω + ( τ ) and θ ( τ ) + 2 nπ intersect for τ ∈ (0 , τ ∗ ), where n is anonnegative integer.For each integer n >
0, denote the j n points of intersection of the curves τ ω + ( τ )and θ ( τ ) + 2 nπ for τ ∈ (0 , τ ∗ ), in increasing order, by τ jn , j = 1 , . . . , j n , i.e., for each n = 0 , , , . . . , τ jn ω + ( τ jn ) = θ ( τ jn ) + 2 nπ, j = 1 , , . . . , j n . Theorem 3.4.
Consider system (5). The characteristic equation (10) has a pair ofpure imaginary eigenvalues, if and only if, τ = τ jn ∈ (0 , τ ∗ ) , a point of intersection ofthe curves τ ω + ( τ ) and θ ( τ ) + 2 nπ , for some integer n > . At all such intersections,the pair of pure imaginary eigenvalues is simple and no other root of (10) is an integer
20 40 60 80 1000246810 τω + ( τ ) θ ( τ ) + 2 πθ ( τ ) τ τω + ( τ ) θ ( τ ) + 6 πθ ( τ ) + 4 πθ ( τ ) + 2 πθ ( τ ) τ Figure 1.
Intersections of θ ( τ ) + 2 nπ and τω + ( τ ) , n = 0 , , . . . , . Values of τ at which the characteristicequation has pure imaginary eigenvalues, and hence candidates for critical values of τ at which there couldbe Hopf bifurcations. In both graphs, at all such intersections, transversality holds, since the slope of thesecurves at these intersections are different. Parameters: m = 1 , r = 1 , K = 1 , Y = 0 . . (LEFT) s = 0 .
02. For n = 0 there are two intersections (i.e. j = 2), at τ and τ , but for n = 1, and hence n >
1, there are nointersections. (RIGHT) s = 0 . j n = 2 , n = 0 , , τ n and τ n , for n = 0 , n = 3, and hence n >
3, there are no intersections. In both (LEFT) and (RIGHT) , E + is asymptotically stable for τ ∈ [0 , τ ) ∪ ( τ , τ c ) and unstable for τ ∈ ( τ , τ ). multiple of iω + ( τ jn ) . If in addition, dd τ ( τ ω + ( τ )) (cid:12)(cid:12)(cid:12) τ = τ jn = dd τ θ ( τ ) (cid:12)(cid:12)(cid:12) τ = τ jn , the transversalitycondition for Hopf bifurcation, dd τ Re( λ ( τ )) (cid:12)(cid:12)(cid:12) τ = τ jn = 0 , holds. The proof is given in A.6.
Remark 5.
If the slope of the curve θ ( τ ) + 2 nπ is less than the slope of τ ω + ( τ ) atan intersection point τ jn , then a pair of complex roots of (10) crosses the imaginaryaxis from left to right as τ increases through τ jn . On the other hand, if the slope ofthe curve θ ( τ ) + 2 nπ is greater than the slope of τ ω + ( τ ) at an intersection point τ jn ,then a pair of complex roots of (10) crosses the imaginary axis from right to left as τ increases through τ jn . Corollary 3.5.
Consider system (5). Assume that τ ∈ [0 , τ ∗ ] and that there exists N > such that (2 N + 1) π max τ ∈ [0 ,τ ∗ ] τ ω + ( τ ) N + 1) π .(1) For n N , θ ( τ ) + 2 nπ and τ ω + ( τ ) have at least two intersections in (0 , τ ∗ ) .(2) For n > N + 1 , θ ( τ ) + 2 nπ and τ ω + ( τ ) do not intersect in (0 , τ ∗ ) .(3) If θ ( τ ) + 2 nπ and τ ω + ( τ ) intersect for any n > , then τ is the smallest and τ j the largest value of τ for which (10) has a pair of pure imaginary eigenvalues.(4) The coexistence equilibrium E + is locally asymptotically stable for τ ∈ [0 , τ ) ∪ ( τ j , τ c ) . The proof is given in A.7.Our results differ from those in Gourley and Kuang [12], since we give explicitformulas for solutions of (13) and define θ ( τ ) explicitly. These explicit formulas playan important role in analysis and make numerical simulations more straightforward.Although ω − ( τ ) in equation (14) is negative in this model and hence its square rootis not real, in other models equation (13) can have two positive real solutions (see [7],where both ω + ( τ ) and ω − ( τ ) are positive). In that case, double Hopf bifurcations arepossible.In the next section we will demonstrate numerically that in the example shownin Figure 1, E + first loses its stability through a supercritical Hopf bifurcation as τ τ and then restabilizes as a result of a second supercritical Hopfbifurcation as τ increases through τ . We will also show that between these two valuesof τ there is a sequence of bifurcations resulting in interesting dynamics, including astrange attractor.
4. An example demonstrating complex dynamics
In this section, unless specified otherwise, we select the following values for the pa-rameters in model (2): m = 1 , r = 1 , K = 1 , Y = 0 . , s = 0 . , (21)and consider τ as a bifurcation parameter. Since, with this selection of parameters, themodel is already in the form of the scaled version of the model (5), it is not necessaryto apply the scaling given by (4). We first use this example to illustrate the analyticresults given in section 3 where we provided necessary and sufficient conditions for asimple pair of pure imaginary eigenvalues of the characteristic equation to occur as τ varies. We then provide bifurcation diagrams with τ as the bifurcation parameter,simulations including time series and time delay embeddings for various values of τ , anda return map at a value of τ at which there is a chaotic attractor, in order to illustratethe wide variety of dynamics displayed by the model, even in the case when there areonly two Hopf bifurcations. This includes two supercritical Hopf bifurcations, saddle-node bifurcations of limit cycles and sequences of period doublings that appear to leadto chaotic dynamics with a strange attractor reminiscent of the strange attractor forthe well-known Mackey-Glass equation [20, 21] dxdt = β x ( t − τ )1 + ( x ( t − τ )) n − γx ( t ) . For the parameters given by (21), the model has three equilibria: E = (0 , E =(1 , E + = ( x + ( τ ) , y + ( τ )) = (0 . e . τ , − . e . τ ), givenby (7). The components of E + are both positive, if, and only if, τ ∈ [0 , τ c ), where τ c ≈
170 (see (8)), and by part 3(b) of Theorem 3.1 the model is uniformly persistentfor τ ∈ [0 , E is always a saddle, and hence unstable. E is globally asymptotically stablefor τ > τ c , and unstable for τ ∈ [0 , τ c ). For τ = 0, E + is asymptotically stable,and by Lemma 3.2, can only lose stability by means of a Hopf bifurcation. In order todetermine th:e critical values of τ at which there is a pair of pure imaginary eigenvalues: λ ( τ ) = ± iω ( τ ), we consider the interval (0 , ω ( τ ) is positive, is bounded above by τ ∗ ≈
115 (defined in (15)).With the parameters given in (21), by (16) the function ω + ( τ ) = vuut − . ˙5 10 − ( e . τ ) + q .
235 10 − ( e . τ ) + . − e . τ − . e . τ + 4210nd recall, by (19) θ ( τ ) = arccos(cos( τ ω + ( τ ))) . The intersections of the functions τ ω + ( τ ) and θ ( τ ) + 2 nπ, n a nonnegative integer, inthe interval (0 , τ ∗ ), give the critical values of τ for which there is a pair of pure imag-inary eigenvalues. There are only two intersections as can be seen in Figure 1 (LEFT) .Since, π < max τ ∈ [0 ,τ ∗ ] τ ω ( τ ) < π , by Corollary 3.5, we are guaranteed at least twovalues of τ > (LEFT )), there are precisely two such values, τ and τ . UsingMaple [22], we found that τ ≈ .
917 and τ ≈ . θ ( τ ) and τ ω + ( τ ) are different at these intersection points, since thesecurves cross transversally (see Figure 1 (LEFT) ), and hence the transversality requiredfor Hopf bifurcation holds at each root. By part 4 of Corollary 3.5 and Remark 5, thestability of E + changes from asymptotically stable to unstable as τ increases through τ and from unstable to asymptotically stable as τ increases through τ , and is un-stable for τ ∈ ( τ , τ ). But recall, even though E + is unstable here, by part 3(b) ofTheorem 3.1, the model is still uniformly persistent in this interval.Thus, we have shown that for the parameters chosen, there are exactly two candi-dates for Hopf bifurcations (see [26], Chapter 6, Theorem 6.1, page 89-90). That bothHopf bifurcations are supercritical (involving first the birth, and then the disappear-ance of orbitally asymptotically stable periodic solutions) will be demonstrated in thenext section. As τ increases through τ = 1 . τ increases, and that they disappear when τ increases throughthe critical value τ = 108 . s , the value of n for which the curves θ ( τ ) + 2 nπ and τ ω + ( τ ) inter-sect can increase. See Figure 1 (RIGHT) for an example with s = 0 .
007 for which thecurves θ ( τ ) + 2 nπ and τ ω + ( τ ) intersect when n = 0 , ,
2. In fact, there are 6 pointsof intersection. We see by Remark 5, that for this example, E + is asymptotically sta-ble until τ = τ , is unstable for τ ∈ ( τ , τ ), and finally becomes stable again for τ ∈ ( τ , τ c ). Again, even though the model is unstable for τ ∈ ( τ , τ ), by part 3(b) ofTheorem 3.1, it is uniformly persistent in this interval, since the model is uniformlypersistent whenever both components of E + are positive, i.e. for τ ∈ [0 , τ c ). The computations and figures in this section were done using Maple [22], MATLAB[23], and XPPAUT [6]In Figure 2 (TOP) , for each value of τ ∈ [0 , x ( t ) = y ( t ) = 0 . t ∈ [ − τ, y ( t ) coordinate on the attractor. Because we areinterested in period doubling bifurcations, we then eliminate certain local maximaand minima that are due to kinks in the solutions (see Figure 3) rather than actualbifurcations, to obtain the graph in Figure 2 (BOTTOM) ). Our solutions have kinks forvalues of τ ∈ [55 , τ ≈ .
917 and τ ≈ .
20 40 60 80 100 12000.511.522.5 τ y Figure 2.
Orbit diagrams. Initial data was taken to be x ( t ) = y ( t ) = 0 . t ∈ [ − τ, (TOP) All local maxima and minima for the y ( t ) coordinate of the attractoras τ varies, including kinks. (BOTTOM) Diagram including local maxes and mins for the y ( t ) coordinate as τ varies, but with kinks eliminated. There are two saddle-node of limit cycle bifurcations. They occur for τ approximately equal to 76 and 82, where the curves in the orbit diagrams stop abruptly and there appear tobe vertical dots. For τ between these values, there is is bistability. Two orbitally asymptotically stable periodicorbits (with their maximum and minimum amplitudes shown) and an unstable periodic orbit with amplitudesbetween them (not shown). The two stable periodic orbits were found by producing this part of the orbitdiagram varying τ forward and then varying it backwards but startng at the last point of the attractor for theprevious value of τ . τ=70 kink kink Figure 3.
The time series for y ( t ) when τ = 70, depicting kinks. There are two local maxima and two localminima over each period as shown in Figure 2 (TOP), but only one local maxima and one local minima inFigure 2 (BOTTOM) in which kinks have been removed.
75 80 85 90 95 10000.511.5
Figure 4.
Zoom-in of orbit diagram shown in Figure 2 for τ ∈ [75 , τ ∈ [80 , τ ≈ . τ and another saddle-node bifurcation oflimit cycles for a value of τ smaller than τ = 81. For values of τ between these twosaddle-node bifurcations, there is bistability. There are two orbitally asymptoticallystable period orbits. An example of two such orbits is given in Figure 5, where τ = 81. τ=81 y(t- τ) Figure 5.
Time delay embedding of two orbitally asymptotically stable periodic orbits demonstrating bista-bility for τ = 81. The one with larger amplitude (dashed) has initial data x ( t ) = y ( t ) = 0 .
1, for t ∈ [ − τ, x ( t ) = y ( t ) = 0 .
1, for t ∈ [ − τ,
0) and x (0) = 0 . , y (0) = 0 .
83 and has period approximately 273 . Figure 4 suggests that there are sequences of period doubling bifurcations, oneinitiating from the left at τ ≈
83, 86, 86 . . . . , and one initiating from the right at τ ≈ .
3, 93 .
2, 92 .
2, and τ between 92 and 91 .
85. To demonstrate these sequences,time series ( y ( t ) versus t ) and time delay embeddings ( y ( t ) versus y ( t − τ )) at valuesof τ between these bifurcations are shown in Figures 6 and 7.Figure 4 also suggests that between these sequences of period doubling bifurcationsthere is a window of values of τ at which there are periodic attractors that do nothave a period that results from a bifurcation with period approximately equal to 2 n for some integer n , as well as chaotic dynamics. An example of the former is illustratedin Figure 8. For τ = 90 .
7, the time-series embedding of a periodic orbit with periodapproximately equal to 1800 time steps involving six loops (2 times 3) is shown.The time series and the time delay embedding of a chaotic attractor for τ = 90is shown in Figure 9. The projection of this attractor into ( x ( t ) , y ( t ))-space is alsoshown in Figure 10. This strange attractor resembles the chaotic attractor of the well-known Mackey-Glass equation ([21], Figure 2). The return map shown in Figure 11,for τ = 90, also resembles the return map for the Mackey-Glass equation ([21], Figure14) in the case of chaotic dynamics. Sensitivity to initial data is a hallmark of chaoticdynamics. Figure 12 (RIGHT) demonstrates that there is sensitivity to initial data inthe case of the solution for τ = 90 that converges to the strange attractor, shown inFigure 9. To show that this is not just a numerical artifact, in Figure 12 (LEFT) we showthat, as expected, there is no sensitivity for the solution for τ = 92 that converges tothe periodic solution shown in Figure 7.This example demonstrates that including delay in a simple predator-prey modelthat always has a globally asymptotically stable equilibrium point in the absence ofdelay, cannot only destabilize a globally asymptotically stable equilibrium point, butcan even result in the birth of a strange attractor.14 τ=82 τ=82 y(t- τ) τ=85 τ=85 y(t- τ) τ=86.3 τ=86.3 y(t- τ) τ=86.8 y(t- τ) τ=86.8 Figure 6. (LEFT)
Time series starting from the initial data x ( t ) = y ( t ) = 0 . , t ∈ [ − τ,
0] indicating howquickly the orbit gets close to the periodic attractor and (RIGHT) time delay embeddings of the periodicattractors, demonstrating the sequence of period doubling bifurcations initiating from the left at τ ≈ , .
6. Values of τ selected between these bifurcations: τ = 82 , , .
3, and 86 .
8, with periods of theperiodic attractor approximately equal to: 340 . , . , .
5, and 2298 .
3, respectively, are shown. τ=91.95 τ=91.95 y(t- τ) τ=92 y(t- τ ) y(t) τ=92 τ=93 τ=93 y(t- τ) τ=96 τ=96 y(t- τ) τ=100 τ=100 y(t- τ) Figure 7. (LEFT)
Time series starting from the initial data x ( t ) = y ( t ) = 0 . , t ∈ [ − τ,
0] indicating howquickly the orbit gets close to the periodic attractor and (RIGHT) time delay embeddings of the periodicattractors, demonstrating the sequence of period halfing bifurcations initiating from the left for values of τ between 91 . τ ≈ . , . , .
3. Graphs shown are for values of τ between these bifurcations: τ = 91 . , , ,
96, and 100, with periods approximately equal to: 4794 . , . , . , . , and295 .
3, respectively. τ =90.7 y(t- τ ) Figure 8.
Time delay embedding starting at initial data x ( t ) = y ( t ) = 0 . , t ∈ [ τ,
0] for τ = 90 . = 2 n ) loops for some integer n . τ=90 τ=90 y(t- τ) Figure 9. (LEFT)
Time series for τ = 90 starting from the initial data x ( t ) = y ( t ) = 0 . , t ∈ [ − τ, (RIGHT) Time delay embedding of the strange attractor for τ = 90. Only the portion of the orbit from t = 240 , − ,
000 is shown. τ=90
Figure 10.
The strange attractor, for τ = 90, shown in Figure 9 in ( x, y )-space. Only the portion of theorbit from t = 240 ,
000 to 260 ,
000 is shown. min − y m i n τ = 90 Figure 11.
The return map for τ = 90, computing the minimum value of y ( t ) as a function of the precedingminimum value of y ( t ) for y ( t ) < . Figure 12.
Time series (LEFT) for the solution that converges to a periodic attractor when τ = 92, and (RIGHT) for the solution that converges to a strange attractor when τ = 90, demonstrating that there is nosensitivity to initial data in the former case, but that there is sensitivity in the latter case. Initial data usedfor the solid curves: x ( t ) = y ( t ) = 0 . t ∈ [ − τ, x ( t ) = 0 .
11 and y ( t ) = 0 . t ∈ [ τ, . Discussion and Conclusions We investigated the effect of the time required for predators to process their prey onthe possible dynamics predicted by a mathematical model of predator-prey interaction.We incorporated a discrete delay to model this process in one of the simplest classicalpredator-prey models, one that only allows convergence to an equilibrium when thisdelay is ignored. We showed that including the delay results in a model with muchricher dynamics. By choosing one of the simplest models when delay is ignored, onethat predicts that no oscillatory behaviour is possible, the effect of the delay on thedynamics is emphasized.This model can also be interpreted as a model of a stage-structured population withthe delay modelling the maturation time of the juveniles (see Gourley and Kuang [12]).In the model we considered, the prey are assumed to grow logistically in the absenceof the predator. The interaction of the predator and prey is described using a linearresponse function, often referred to as mass action or Holling type I. It is well-knownthat when delay is ignored this is one of the simplest predator-prey models for whichall solutions converge to a globally asymptotically stable equilibrium point for allchoices of the parameters. Therefore, any resulting non-equilibrium dynamics wouldthen be solely attributable to the introduction of the delay in the growth term of thepredator. We not only found non-trivial periodic solutions, but also bistability, andchaotic dynamics. It is then likely that there is similar rich dynamics in most predator-prey models with any reasonable response function, when such a delay is incorporated,for some selection of the parameters, including the form most used by ecologists, theHolling type II form. This form given mathematically by f ( x ) = mx/ (1 + bx )), canbe considered a generalization of the Holling type I form, obtained by simply addingan extra parameter, b . However, we feel that demonstrating that this wide range ofdynamics is possible for even for one of the simplest models gives more compellingevidence that delay should not be ignored when making policy decisions.Understanding how changes in average temperature might affect survivability ofendangered populations or result in invasions by undesirable populations is important.Since temperature can affect how quickly predators process the prey that they capture,based on our results there might be important implications for populations in thewild. In most predator populations, the processing time, τ , is faster when it is warmerand slower when it is colder. Our results may help us understand how a change inaverage temperature might influence the dynamics of particular predator-prey systemsof interest.Our study suggests that we need to be careful when measuring population sizesin the wild and predicting the general health of the population based upon whetherthe population seems to be increasing or decreasing. Short term indications that apopulation size is changing in a system with oscillatory dynamics may be misleadingand predicting future population size may be impossible without more information.It would be necessary to have some idea of the period of the intrinsic oscillations ifthe population is suspected to vary periodically. If the dynamics are suspected to bechaotic, this may be even more complicated, due to sensitivity to initial data.If a predator-prey system has the potential to have chaotic dynamics, based on ourresults, is there anything that is predictable? Can such analyses suggest how to preventextinctions or invasions due to a change in average temperature that could result in achange in the processing time of the prey by the predator? We give some observationsbased on the predictions of our model and the example we considered in Section 4.However, more work would need to be done to determine whether these predictions19re consistent for more realistic models, and if so, long term observations would haveto be made by ecologists to determine if they are relevant for populations in the wild.First, at the one extreme, if the processing time is too long, the predator populationwould not be expected to survive, since it is obvious that if the processing time is longerthan the life span of the predator, the predator population has no chance to avoidextinction. If there was no such threshold for extinction predicted by the model, themodel should be abandoned. It is therefore very important to use the term e − sτ y ( t − τ )and not just y ( t − τ ) in order to account for the predators that do not survive longenough to affect growth of the population, in the equation describing the growth ofthe predator in model (2). Due to this term, in our model there is such a threshold, τ c .From the orbit diagram (see Figure 2), for the example considered in Section 4,we summarize some observations. The dynamics for both populations are oscillatoryfor a wide range of processing times, and non-oscillatory for only relatively (very)short or relatively long processing times (i.e. before the first Hopf bifurcation at τ ≈ . τ ≈ τ at the final Hopf bifurcation of the coexistence equilibrium,the population is no longer oscillatory. The size of the predator population decreasesrelatively slowly as the processing time increases further. Although the size of thepredator population gets smaller as the processing time increases, it does not get muchsmaller, until the processing time gets close to the threshold for extinction τ c . If wehad shown the diagram extended to the threshold τ c = 170, one would see that as theprocessing time gets close to τ c , the size of the predator population suddenly decreasesrelatively quickly to zero. So our model predicts that for a predator population withfairly long processing times to begin with, cooling of the environment could be expectedto be detrimental with respect to the survivability of the predator population. Thus,this effect on the size of the predator population might be minor if the delay is closeto its Hopf bifurcation value, but could be drastic if it is close to the extinction value.Our example also suggests that the predator population may not be oscillatorywhen it becomes endangered and hence close to extinction, i.e., for excessively longprocessing times. It would be interesting to investigate if this also holds for the model,with Holling type I response functions replaced by Holling type II. In the model withHolling type II response functions, if the carrying capacity of the environment for theprey is not affected by the cooling and it is relatively high, then the model predictsthat the predator population could still be oscillatory, even if the processing timefor the predator is ignored. However, cooling might also be expected to reduce thecarrying capacity for the prey, moving the parameters to a range where that modelwould also predict non-oscillatory dynamics (the paradox of enrichment [25]). So onceagain, perhaps alarm bells should be sounded when there is cooling and the predatorpopulation has been non-oscillatory and appears to decline rapidly as the averageyearly temperature declines.On the other hand, at the other extreme, our model predicts a Hopf bifurcationat a relatively small value of the delay, resulting in the birth of a family of periodicorbits with amplitude increasing very quickly as the delay increases with both the preyand predator populations spending time very close to zero. This remains the case forsmall, but intermediate values of the delay, (before the two saddle node of limit cyclesbifurcations near τ = 80). For this range of τ , these populations are therefore verysusceptible to stochastic extinctions. If the processing time was originally very small(below the critical value for the first Hopf bifurcation) and cooling made it longer,again a stochastic extinction might be likely. Similarly, if τ was close to the first of20he two saddle node of limit cycles bifurcations near τ = 80), warming of the averagetemperature could result in a stochastic extinction of one or both of the populations.Since the predator cannot survive without the prey in our model, even if it was theprey population that experienced the stochastic extinction, the predator populationwould eventually die out as well. Between the smallest value of τ at which there is aperiod doubling bifurcation, and the value of τ at the largest Hopf bifurcation, coolingwould probably be advantageous to the predator population size.In summary, it seems that whether cooling is beneficial or detrimental to the sizeof the predator population depends on where the delay is on the bifurcation diagram,and it is very likely that this is very difficult to determine. As well, if values of τ were to lie in the chaotic region, then changes in the environment that changed thevalue of τ to obtain regular oscillatory dynamics may or may not be preferable. Aswell, as can be seen by the various attractors shown in Figures 6-10, and the orbitdiagrams in Figure 4, certain properties of the system in the chaotic region such asthe maximum and minimum values of the predator size were fairly insensitive to thechange in the delay. This is only a toy model. However, it suggests that ignoring delayin a model can result in incorrect predictions. Here, delay could change the dynamicsfrom convergence to a globally asymptotically stable equilibrium to wild oscillations.It is most likely that in nature, it would not be possible to distinguish from data, apriori, if the dynamics were chaotic or periodic, but more importantly it would notnecessarily be predictable what the effect of an increase or a decrease in the delaywould be. Hence, our analysis indicates that one should be extra cautious if trying tomanipulate the delay to achieve a certain result based on model predictions.Finally it is worth pointing out that the resulting strange attractor in the predator-prey model studied here bares such a close resemblance to the Mackey-Glass attractor,a model involving a single delay differential equation to model a simple feedback systemfor respiratory control or hematopoietic diseases. Understanding whether there is adeeper significiance to this relationship may give us a better understanding of theclass of possible strange attractors and warrants further investigation. Funding
The Research of Gail S. K. Wolkowicz was partially supported by Natural Sciencesand Engineering (NSERC) Discovery Grant
References [1]
J. Arino, L. Wang, and G. S. K. Wolkowicz , An alternative formulation for adelayed logistic equation , J. Theoret. Biol., 241(1) (2006), pp. 109–119.[2]
A. D. Bazykin , Nonlinear Dynamics of Interacting Populations , vol. 11 of A, WorldScientific Series on Nonlinear Science, Singapore, 1998.[3]
E. Beretta and Y. Kuang , Geometric stability switch criteria in delay differentialsystems with delay dependent parameters , SIAM J. Math. Anal., 33 (2002), pp. 1144–1165.[4]
K. L. Cooke, R. H. Elderkin, and W. Huang , Predator-prey interactions with delaysdue to juvenile maturation , SIAM J. Appl. Math., 66 (2006), pp. 1050–1079.[5]
R. Driver , Existence and stability of solutions of a delay-differential system , Arch. Ra-tion. Mech. Anal., 10 (1962), pp. 401–426. B. Ermentrout , Simulating, Analyzing, and Animating Dynamical Systems: A Guideto XPPAUT for Researchers and Students , SIAM, 2002.[7]
G. Fan and G. S. K. Wolkowicz , A predator-prey model in the chemostat with timedelay , Int. J. Differ. Equ., (2010), pp. Art. ID 287969, 41.[8]
J. E. Forde , Delay differential equation models in mathematical biology , PhD thesis,University of Michigan, 2005.[9]
H. I. Freedman , Deterministic Mathematical Models in Population Ecology , MarcelDekker, New York, 1980.[10]
F. R. Gantmacher , Applications of the Theory of Matrices , Trans. J. L. Brenner et al.,New York: Interscience, 1959.[11]
A.-M. Ginoux, B. Rossetto, and J.-L. Jamet , Chaos in a three-dimensionalVolterra-Gause model of predator-prey type , Internat. J. of Bifur. Chaos Appl. Sci. Engrg.,15 (2005), pp. 1689–1708.[12]
S. A. Gourley and Y. Kuang , A stage structured predator-prey model and its depen-dence on maturation delay and death rate , J. Math. Biol., 49 (2004), pp. 188–200.[13]
J. K. Hale and V. L. S. M. , Introduction to Functional-Differential Equations , vol. 99of Applied Mathematical Sciences, Springer-Verlag, New York, 1993.[14]
M. A. Haque , A predator-prey model with discrete time delay considering different growthfunction of prey , Adv. Appl. Math. Biosci., 2 (2011), pp. 1–16.[15]
A. Hastings and T. Powell , Chaos in a three-species food chain , Ecology, 72 (1991),pp. pp. 896–903.[16]
W. M. Hirsch, H. Hanisch, and J.-P. Gabriel , Differential equation models of someparasitic infections: methods for the study of asymptotic behavior , Comm. Pure Appl.Math., 38 (1985), pp. 733–753.[17]
Y. Kuang , Delay Differential Equations with Applications in Population Dynamics ,vol. 191 of Mathematics in Science and Engineering, Academic Press Inc., Boston, MA,1993.[18]
V. Lakshmikantham and S. Leela , Differential and Integral Inequalities: Theory andApplications , vol. 55-I, Academic Press, New York, 1969.[19]
M. Y. Li, X. Lin, and H. Wang , Global Hopf branches and multiple limit cycles in adelayed Lotka-Volterra predator-prey model , Discrete Contin. Dyn. Syst. Ser. B, 19 (2014),pp. 747–760.[20]
M. C. Mackey and L. Glass , Oscillations and chaos in physiological control , Science,197 (1977), pp. 287–289.[21] ,
Mackey-Glass equation , Scholarpedia, 4 (2009), p. 6908.[22]
MAPLE , Maplesoft, a division of Waterloo Maple Inc. , Waterloo, Ontario, 2017.[23]
MATLAB , version 9.5.0 (R2018b) , The MathWorks Inc., Natick, Massachusetts, 2018.[24] A. Morozov, S. Petrovskii, and B.-L. Li , Bifurcations and chaos in a predator-preysystem with the allee effect , P Roy. Soc. B-Biol. Sci., 271 (2004), pp. 1407–1414.[25]
M. L. Rosenzweig , Paradox of enrichment: Destabilization of exploitation ecosystemsin ecological time , Science, 171 (1971), pp. 385–387.[26]
H. L. Smith , An Introduction to Delay Differential Equations with Applications to theLife Sciences , vol. 57 of Texts in Applied Mathematics, Springer, New York, 2011.[27]
H. L. Smith and H. R. Thieme , Dynamical Systems and Population Persistence ,vol. 118 of Graduate Studies in Mathematics, American Mathematical Society, Provi-dence, RI, 2011.[28]
J. Wang and W. Jiang , Bifurcation and chaos of a delayed predator-prey model withdormancy of predators , Nonlinear Dynam., 69 (2012), pp. 1541–1558. ppendix A. Proofs A.1. Proof of Proposition 2.1
Proof.
Parts 1 & 2. Assume that ( φ ( t ) , ψ ( t )) ∈ X . Since the right hand side of (5) isLipschitz, there exist h > t ∈ [0 , h ].We will show that solutions are bounded and hence do not blow up in finite time.Therefore, by the standard results for existence and uniqueness of solutions for delaydifferential equations in Driver [5], it will follow that solutions exist and are uniquefor all t ≥ φ (0) >
0, and the face where x ( t ) = 0 is invariant, that x ( t ) > t > y ( t ) > t >
0, follows since y ′ ( t ) > − sy ( t ) for all t ∈ [0 , τ ] and so y ( t ) > t ∈ [ − τ, τ ]. Arguing inductively on [ nτ, ( n + 1) τ ] , n = 1 , , . . . , it follows that y ( t ) > t >
0, where the solution exist.To show that solutions are bounded above, consider the first equation of (5)˙ x ( t ) = x ( t )(1 − x ( t )) − y ( t ) x ( t ) x ( t )(1 − x ( t )) . It is well-known that for the logistic equation ˙ z ( t ) = z ( t ) (1 − z ( t )), given any ǫ > T >
0, such that | z ( t ) | < ǫ for all t > T . Using a comparison principle(e.g. [18] Theorem 1.4.1, page 15) x ( t ) z ( t ), and so 0 x ( t ) < ǫ for all t > T .To prove that y ( t ) is bounded above, define w ( t ) = Y e − sτ x ( t − τ ) + y ( t ) . (A1)Then ˙ w ( t ) = Y e − sτ dx ( t − τ ) dt + dy ( t ) dt , = − sy ( t ) + Y e − sτ x ( t − τ ) (1 − x ( t − τ )) , = − sw ( t ) + Y e − sτ x ( t − τ ) ( s + 1 − x ( t − τ )) , − sw ( t ) + 14 Y e − sτ ( s + 1) , since (cid:0) x ( t − τ ) − s +12 (cid:1) > x ( t − τ ) ( s + 1 − x ( t − τ )) ( s +1) . Therefore,since z ( t ) = z (0) e − st + s Y e − sτ ( s + 1) (1 − e − st ) is the solution of the initial valueproblem ˙ z ( t ) = − sz ( t ) + 14 Y e − sτ ( s + 1) , z (0) = w (0) > , using a comparison principle it follows that w ( t ) z ( t ) for all t >
0. Consequently, by(A1), y ( t ) w ( t ) w (0) e − st + s Y e − sτ ( s + 1) (1 − e − st ). Therefore, y ( t ) is bounded.Part 3. Assume that ( φ ( t ) , ψ ( t )) ∈ X . Since φ (0) >
0, and the face where x ( t ) = 0is invariant, it follows that x ( t ) > t >
0. Since there exists θ ∈ [ − τ,
0] such that φ ( θ ) ψ ( θ ) >
0, take T = τ + θ , and note that 0 T τ . Since y ( T ) >
0, either y ( T ) > y ( T ) = 0 and by (5), y ′ ( T ) = − sy ( T ) + Y e − sτ y ( θ ) x ( θ ) = 0 + Y e − sτ ψ ( θ ) φ ( θ ) > ǫ > y ( t ) > t ∈ ( T, T + ǫ ]. But thisimplies that y ′ ( t ) > − sy ( t ) for all t ∈ [ − τ, T + ǫ + τ ]. Therefore, y ( t ) > ∈ [ T + ǫ, T + ǫ + τ ]. By repeating this argument, it follows that y ( t ) > t > T + ǫ for any ǫ > ✷ A.2. Proof of Theorem 3.1
Proof.
Part 1. Evaluating (9) at E gives P ( λ ) | E = ( λ + s )( λ −
1) = 0 , which hastwo real roots λ = − s and λ = 1. Therefore, E is a saddle.Part 2.(a) Evaluating (9) at E gives P ( λ ) | E = ( λ + 1)( λ + s − Y e − ( s + λ ) τ ) = 0 . One of the roots is λ = −
1. The other roots satisfy g ( λ, τ ) := ( λ + s ) e ( λ + s ) τ = Y. Forany fixed 0 τ < τ c , there is a real root λ ( τ ) >
0, such that g ( λ ( τ ) , τ ) = Y , since g (0 , τ ) < Y and g ( λ, τ ) → ∞ as λ → ∞ . Hence, E is unstable.Part 2.(b) Next assume that τ > τ c . We prove that E is globally asymptoticallystable. Since se sτ Y > ǫ := (cid:0) sY e sτ − (cid:1) >
0, and so by Proposition 2.1, there existsa
T > x ( t ) < ǫ for all t > T . Therefore, for all sufficiently large t , Y e − sτ x ( t − τ ) < Y e − sτ (1 + ǫ ) = Y e − sτ (cid:0) (cid:0) sY e sτ − (cid:1)(cid:1) = Y e − sτ + s < s, and so the second equation of (5) can be written ˙ y ( t ) = − sy ( t ) + b ( t ) y ( t − τ ) , where b ( t ) := Y e − sτ x ( t − τ ) < s . Choosing α = s/ b ( t ) < s = 4( s − α ) α = 4( s − s/ s/ y ( t ) → t → ∞ .Hence, for any ǫ >
0, there exists T such that 0 < y ( t ) < ǫ for t > T . From thefirst equation of (5), for any 0 < ǫ < x ( t ) (1 − x ( t ) − ǫ ) < ˙ x ( t ) < x ( t ) (1 − x ( t )) . Note that all solutions of ˙ z ( t ) = z ( t )(1 − z ( t )) converge to z ( t ) = 1 and all solutions of˙ z ( t ) = z ( t )(1 − z ( t ) − ǫ ) converge to z ( t ) = 1 − ǫ as t → ∞ . By a standard comparisonprinciple, for any solution x ( t ) of (5), x ( t ) → t → ∞ . Therefore, E is globallyasymptotically stable.Part 3. Existence of E + follows from (7).Part 3.(a) When τ = 0, model (5) reduces to a well studied model involving onlyordinary differential equations. That E + is globally asymptotically stable in this caseis well known and can easily be proved using phase plane analysis and the Bendixson-negative criterion to rule out periodic solutions.Part 3.(b) The proof is similar to the approach used in Chapter 5.7 of Smith andThieme [27]. The predator reproduction number at a constant level of x is given by R ( x ) = Y e − sτ s x. Then, since E + exists, R (1) > R ( x + ) = 1 . (A2)Let, x ∞ ≡ lim sup t →∞ x ( t ) and y ∞ ≡ lim sup t →∞ y ( t ). Then, we claim x ∞ ≥ se sτ Y . (A3)Suppose not, i.e., that x ∞ is smaller than this minimum. Applying the fluctuationlemma (see Hirsch et. al. [16] or Smith and Thieme [27]) to the y equation in (5),there exists a monotone increasing sequence of times { t n } → ∞ as n → ∞ suchthat ˙ y ( t n ) = 0 and y ( t n ) → y ∞ as n → ∞ . Therefore, 0 = ˙ y ( t n ) = − sy ( t n ) + Y e − sτ x ( t n − τ ) y ( t n − τ ). Letting n → ∞ , we obtain 0 ≤ ( − s + Y e − sτ x ∞ ) y ∞ , sincelim sup n →∞ x ( t n − τ ) y ( t n − τ ) ≤ x ∞ y ∞ . The term in the brackets in this inequality isnegative, since we are assuming (A3) does not hold. It follows that y ∞ = 0. But then,applying the fluctuation lemma to the x equation in (5), it follows that x ∞ = 1 > se sτ Y ,contradicting our assumption that (A3) is not satisfied. Hence, (A3) holds.Let Φ : R + × X → X denote the semiflow generated by model (5). Consider the24unction ρ : X → R + , defined for ( φ ( t ) , ψ ( t )) ∈ X by ρ (( φ, ψ )) = φ (0) so that ρ (Φ( t, ( φ, ψ ))) = x ( t ). By part 3 of Proposition 2.1, ρ (( φ, ψ ))) > ρ (Φ( t, ( φ, ψ ))) > t ≥
0. Therefore, the semi-flow is uniformly weakly ρ -persistent. By a similar argument to that given in Theorem 5.29 of [27] together withProposition 2.1, Φ is a continuous semiflow that has a compact attractor of boundedsets. Therefore, φ is uniformly ρ -persistent by Theorem 5.2 of [27], and hence thereexists ǫ > x ∞ ≡ lim inf t →∞ x ( t ) ≥ ǫ , for all solutions with x (0) > φ, ψ ) ∈ X , there exists ǫ > t →∞ y ( t ) >ǫ , we use Laplace transforms to show that if ( φ, ψ ) ∈ X , then R ( x ∞ ) ≤
1. We denotethe Laplace transform of a function f ( t ) as b f ( λ ) = R ∞ e − λt f ( t ) dt , and note that theLaplace transform of y ( t ) exists for all λ >
0, since y ( t ) is bounded by Proposition 2.1.Taking the Laplace transform on both sides of the y equation of (5) and simplifyingwe obtain:( λ + s ) b y ( λ ) = y (0) + Y e − ( s + λ ) τ c xy ( λ ) + Y e − ( s + λ ) τ Z − τ e − λτ x ( t ) y ( t ) dt. There exists δ > x ( t ) ≥ ( x ∞ − δ ) for all t ≥
0, and so c xy ( λ ) ≥ ( x ∞ − δ ) b y ( λ ) . ( λ + s ) b y ( λ ) ≥ Y e − ( s + λ ) τ ( x ∞ − δ ) b y ( λ ) . By part 3 of Proposition 2.1, y ( t ) > b y ( λ ) > b y ( λ ) to obtain( λ + s ) ≥ Y e − ( s + λ ) τ ( x ∞ − δ ) . After a shift of time, if necessary, we can take the limit as δ, λ → + to obtain s ≥ Y e − sτ x ∞ . This is equivalent to R ( x ∞ ) ≤ . (A4)Define ρ : X → R + by ρ (( φ, ψ )) = min { φ (0) , ψ (0) } .Then, ρ (Φ( t, ( φ, ψ ))) = min { x ( t ) , y ( t ) } . Suppose ( φ, ψ ) ∈ X and y ( t ) is not uniformlypersistent. By part 3 of Proposition 2.1, shifting time if necessary, there is no lossof generality if we assume that ψ (0) = y (0) >
0. Therefore, ρ (Φ( t, ( φ, ψ ))) > t ≥
0. Since by Proposition 2.1, Φ has a compact attractor of bounded sets, byTheorem 5.2 in [27] we need only show that Φ is uniformly weakly ρ -persistent. Supposenot. Recall that we have already shown that x ∞ > ǫ > . Take any ǫ ∈ (0 , ǫ ). Then,there exists a solution with x (0) > y (0) > y ∞ < ǫ. Apply thefluctuation lemma to the x equation of (5). Therefore, 0 > x ∞ (1 − x ∞ − ǫ ) . Taking ǫ > R (1) > R ( x ) is increasing, it follows that R ( x ∞ ) >
1, contradicting (A4). ✷ .3. Proof of Lemma 3.2: Roots of the characteristic equation cannotbifurcate in from infinity Proof.
In Kuang [17] (Theorem 1.4, Chapter 3, page 66), taking n = 2 and g ( λ, τ ) = p ( τ ) λ + ( qλ + c ( τ )) e − λτ + α ( τ ) , sincelim sup Re λ> , | λ |→∞ | λ − g ( λ, τ ) | = 0 < , no root of (10) with positive real part can enter from infinity as τ increases from 0.Since, when τ = 0 all roots have negative real parts, the result follows. ✷ A.4. Proof of Theorem 3.3
Proof.
Using the quadratic formula to solve for ω in (13), the roots must satisfy ω ± = 12 (cid:18) q − p ( τ ) + 2 α ( τ ) ± p ( q − p ( τ ) + 2 α ( τ )) − α ( τ ) − c ( τ )) (cid:19) . Since q − p ( τ ) + 2 α ( τ ) = s − s (cid:18) e sτ Y (cid:19) + 2 s e sτ Y = − (cid:18) se sτ Y (cid:19) < ,ω − is either complex or negative for any τ >
0, and ω is positive, if, and only if,( α ( τ ) − c ( τ )) = − s (cid:18) se sτ Y − (cid:19) (cid:18) se sτ Y − (cid:19) < . This is the case, if, and only if, se sτ Y < or se sτ Y > . However, E + only exists when x + ( τ ) = se sτ Y <
1. Therefore, a real positive root exists if, and only if, x + ( τ ) = se sτ Y < .This implies that τ < τ ∗ . Hence, for τ ∈ [0 , τ ∗ ), a real positive root ω + ( τ ) exists, andis defined explicitly by (16).If τ = τ ∗ , then ω + ( τ ) = 0, and if τ > τ ∗ , then either ω + ( τ ) is not real or E + doesnot exist. ✷ A.5. Properties of the functions h ( ω, τ ) and h ( ω, τ ) defined in (17) Lemma A.1.
Assume that E + exists and that sY < .(1) If τ ∈ [0 , τ ∗ ] , then h ( ω + ( τ ) , τ ) + h ( ω + ( τ ) , τ ) = 1 .(2) h ( ω + ( τ ∗ ) , τ ∗ ) = 0 . If τ < τ ∗ , then h ( ω + ( τ ) , τ ) > .(3) h ( ω + ( τ ∗ ) , τ ∗ ) = − . If τ < τ ∗ , then − < h ( ω + ( τ ) , τ ) < . Proof.
Part 1. Since h ( ω + ( τ ) , τ ) is equal to the right-hand side of (12a), and h ( ω + ( τ ) , τ ) is equal to the right-hand side of (12b), where ω = ω + ( τ ) satisfies (13),it follows that h ( ω + ( τ ) , τ ) + h ( ω + ( τ ) , τ ) = 1.Part 2. By Theorem 3.3, ω + ( τ ∗ ) = 0, and so h ( ω + ( τ ∗ ) , τ ∗ ) = 0 . Since for t ∈ [0 , τ ∗ ), x + ( τ ) < and by Theorem 3.3, ω + ( τ ) >
0, it follows that the denominator of26 ( τ, ω + ( τ )) is always positive. Hence, h ( ω + ( τ ) , τ ) >
0, if, and only if,0 < s (1 − x + ( τ )) + x + ( τ )(1 − x + ( τ )) − ω ( τ )= s (1 − x + ( τ )) + x + ( τ )(1 − x + ( τ )) − (cid:18) − x ( τ ) + q x ( τ ) + 4 s (3 x + ( τ ) − x + ( τ ) − (cid:19) . This is equivalent to,12 q x ( τ ) + 4 s (3 x + ( τ ) − x + ( τ ) − < s (1 − x + ( τ )) + x + ( τ )(1 − x + ( τ )) . Since x + ( τ ) < , both sides of the above inequality are positive. Squaring, both sidesyields, 14 (cid:0) x ( τ ) + 4 s (1 − x + ( τ ))(1 − x + ( τ )) (cid:1) < s (1 − x + ( τ )) + x ( τ )(1 − x + ( τ )) + 2 sx + ( τ )(1 − x + ( τ )(1 − x + ( τ )) . But this is equivalent to,0 < s x + ( τ )(1 − x + ( τ )) + x ( τ )(1 − x + ( τ ))(1 − x + ( τ ))+ 2 sx + ( τ )(1 − x + ( τ )(1 − x + ( τ )) . This last inequality is satisfied, since by part 2 of Theorem 3.3, x + ( τ ) < , and so h ( ω + ( τ ) , τ ) > τ ∈ [0 , τ ∗ ).Part 3. That h (( ω + ( τ ∗ ) , τ ∗ ) = −
1, follows from a straightforward calculation, aftersubstituting ω + ( τ ∗ ) = 0 and x + ( τ ∗ ) = in (17b). For τ ∈ [0 , τ ∗ ), the result followsimmediately, by parts 1 and 2. ✷ A.6. Proof of Theorem 3.4
Proof.
Assume that ¯ τ ∈ (0 , τ ∗ ). By Remark 3, (10) has a pair of pure imaginaryeigenvalues, if and only if, ¯ τ = τ jn ∈ (0 , τ ∗ ), for some n > , j j n . Since anecessary condition for a root of (10) to exist is that (13) holds, and hence ω (¯ τ ) = ω + (¯ τ ) given by (16), there can be at most one pair of pure imaginary roots for eachsuch τ = ¯ τ , and hence if any such roots exist, they are simple, and no other root of(10) is an integer multiple.In Beretta and Kuang, [3] (Theorem 4.1, equation (4.1) p.1157), it is shown thatsign (cid:18) dd τ Re( λ ( τ )) (cid:19) (cid:12)(cid:12)(cid:12) τ = τ nj = sign (cid:18) dd τ (cid:20) τ ω + ( τ ) − ( θ ( τ ) + 2 nπ ) ω + ( τ ) (cid:21)(cid:19) (cid:12)(cid:12)(cid:12) τ = τ jn . After differentiating, and noting that ω + ( τ jn ) > τ jn ω + ( τ jn ) − ( θ ( τ jn ) +2 nπ ) = 0, it is easy to see that the term on the right has the same sign as27 dd τ [ τ ω + ( τ ) − ( θ ( τ ) + 2 nπ )] (cid:1) (cid:12)(cid:12) τ = τ jn . It follows that transversality holds whenever thegraphs of τ ω + ( τ ) and θ ( τ ) have different slopes at τ jn . ✷ A.7. Proof of Corollary 3.5
Proof.
Parts 1, 2 and 3 follow immediately from Theorem 3.4, since for τ ∈ [0 , τ ∗ ],both curves are continuous, by part 3 of Lemma A.1 and Remark 4, θ ( τ ) ∈ (0 , π ), and τ ω + ( τ ) > τ = 0 and τ = τ ∗ . (See Figure 1 for typical examples.)Next consider local stability of E + when it exists, i.e. τ ∈ [0 , τ c ). By part 3(a) ofTheorem 3.1, E + is globally asymptotically stable when τ = 0. By Lemma 3.2, E + can only change stability for 0 τ < τ c , if a pair of roots cross the imaginary axis. As τ increases from 0, by part 3, the first possible such crossing occurs for τ = τ . Hence, E + is locally asymptotically stable for τ ∈ [0 , τ ). When τ = τ c , E + and E coalesce,and by part 2(b) of Theorem 3.1, E is globally asymptotically stable for τ > τ c . Sinceagain by Lemma 3.2, E + can only change stability for 0 τ < τ c , if a pair of rootscross the imaginary axis, and τ j is the last possible such crossing, E + must be locallyasymptotically stable for τ ∈ ( τ j , τ c ). ✷✷