A nonlinear continuous-time model for a semelparous species

Mathematical Biosciences 297 (2018) 1–11

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

A nonlinear continuous-time model for a semelparous species

T

A. Veprauskas Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504-1010, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Semelparous species Dynamic dichotomy Hopf bifurcations McKendrick model

Periodical semelparous insects such as cicadas and May beetles exhibit synchronization in age classes such that only one age class is present at any point of time. This leads to outbreaks of adults as they all reach maturity around the same time. Discrete-time models of semelparous species have shown that this type of synchronous cycling can occur as a result of greater between-class competition relative to within-class competition. However, relatively few studies have examined continuous-time models of semelparous species. Here we develop a continuous-time model for a semelparous species using a technique called the linear chain trick to convert a nonlinear McKendrick partial differential equation into a finite system of ordinary differential equations. We represent semelparity by a birth function whose age distribution can be made arbitrarily narrow. We show that a Hopf bifurcation may occur in this model as a result of competition between reproducing and non-reproducing classes. This bifurcation leads to stable cycles in which the two classes are out of phase, thus providing a continuous-time counterpart to the synchronous cycles that occur in discrete-time models.

1. Introduction

its boundary are invariant. For the two-stage Leslie semelparous model, consisting of a single juvenile and adult stage, a dynamic dichotomy exists between these two possible dynamics [7,11]. When the extinction equilibrium destabilizes, if both steady states bifurcate forward, then one of these steady states is stable and the other is unstable. Which of the two dynamics is stable depends on the relative levels of competition in the population. If the within-class competition is stronger than the between-class competition, then the equilibria are stable and vice versa. For three dimensions, the stability properties become more complicated due to the increased complexity of the dynamics on the boundary of the positive cone which now has invariant loops that result in dynamically complicated cycle chains [7]. For dimensions greater than three, the dynamics have not been fully determined. While the stability of synchronous cycles has not been establish for higher dimensions, a dichotomy between positive equilibria and the boundary of the positive cone, depending on the within-class and between-class competition, still holds [8]. Relatively few papers have investigated whether equivalent dynamics occur in a continuous-time semelparous model. In a continuous setting, the McKendrick model (also called the Lotka-von Foerster model) is often used to study structured populations [1,18,27]. This linear model is discussed in detail in [6,22]. Gurtin and MacCamy [19] introduced the nonlinear McKendrick model given by the equations

A semelparous species is one in which individual organisms reproduce once in their lifetime and, typically, die soon after. Classic examples of semelparous species include the Pacific salmon and annual plants. Semelparity can also be found in various insects including certain species of cicadas, mayflies, and May beetles; longer lived plants like agave and certain types of bamboo; some species of bony fish, frogs, lizards, and even a few marsupials. Much of the mathematical interest in semelparous species has been motivated by periodical semelparous insects such as cicadas and May beetles. These species exhibit synchronization in age classes such that only one age class is present at a time. As a result, the population dynamics can be described by periodic cycles of length equal to the maturation period and all the adults emerge at the same time. It has been suggested that this synchronization may be a result of strong competition between age classes [4]. Discrete-time Leslie matrix models have been used extensively to study the dynamics of semelparous species [2–4,7–15,23–25,31]. These models have provided mathematical verification that strong betweenclass competition can lead to stable periodic cycles. In particular, semelparous Leslie models exhibit two contrasting dynamics. Either equilibria in which all age classes are present or synchronized cycles in which age classes are separated temporally are stable. Mathematically, these dynamics arise because the projection matrix of a semelparous Leslie model is imprimitive and, as a result, both the positive cone and

E-mail address: [email protected]. https://doi.org/10.1016/j.mbs.2018.01.003 Received 28 June 2017; Received in revised form 5 January 2018; Accepted 8 January 2018 Available online 09 January 2018 0025-5564/ © 2018 Elsevier Inc. All rights reserved.

P (t ) =

∫0

∞

ρ (t , a) da

(1a)

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

∂ρ (t , a) ∂ρ (t , a) = −μ (P (t ), a) ρ (t , a) for t , a > 0 + ∂a ∂t

ρ (t , 0) =

∫0

∞

β (P (t ), a) ρ (t , a) da

ρ (0, a) = ρ0 (a),

positive equilibria that bifurcate from the extinction equilibrium. In Section 4.2, we determine the condition necessary for the positive equilibria to destabilize. In Section 5 we show that, for a semelparous species, if the between-class competition is stronger than the withinclass competition, then it is possible for the positive equilibria to destabilize via a Hopf bifurcation. Finally, in Section 6 we present numerical examples to corroborate these results.

(1b) (1c) (1d)

where P(t) is the total population at time t, ρ(t, a) is the density function defined so that

∫a

a2

2. Model equations

ρ (t , a) da

In order to apply the linear chain trick to the McKendrick model (1), we assume

1

is the number of individuals with ages between a1 and a2 at time t, β(P, a) is the birth function, and µ(P, a) is the death function. The boundary condition (1c), often called the renewal equation, describes the birth process and Eq. (1d) gives the initial age distribution. Since individuals are born with age, a = 0 Eq. (1c), the hyperbolic partial differential equation given by Eq. (1b) assumes that changes in population numbers are only due to death for ages a > 0. In [16,29], a semelparous version of the McKendrick model (1) was studied by assuming a Dirac delta fertility function. Dilão and Lakmeche considered the linear McKendrick model and showed that unbounded oscillations can occur for a semelparous species [16]. Rudkicki and Wierczorek studied a nonlinear model [29]. They assumed that maturation occurs at some fixed age, without loss of generality assume age a = 1, and that the fertility function is given by a Dirac delta function multiplied by a density-dependent birth rate β (P , a) = b (P ) δ (a − 1) . They showed that solutions to this model converge weakly to a periodic function if the initial age distribution has sufficiently small support. While competition induced cycles cannot occur in model (1) under these assumptions, it should be noted that, since the density dependent birth rate is assumed to be a function of the total population P, model (1) does not allow for the distinction between the competitive effects of different stages. The goal of this paper is to develop a nonlinear, continuous-time semelparous model such that the nonlinearity allows for the distinction between the reproducing class and the non-reproducing class. With this model formulation, competition can be categorized as either withinclass or between-class competition, and we can investigate whether or not dynamics similar to the discrete-time semelparous model can occur in a continuous-time model. In particular, we are interested in determining if strong between-class competition can result in stable cycles in which the reproducing and non-reproducing classes are out of sync. Such a dynamic would provide a continuous-time counterpart to the dynamic dichotomy observed in the two-stage semelparous Leslie model. Here we use a technique popularized by MacDonald called the linear chain trick to convert the nonlinear McKendrick model into a finite system of ordinary differential equations [26]. This technique has been used to study the dynamics of continuous-time nonlinear McKendrick models for iteroparous species in [19] with generalizations given in [21]. We assume that the birth function is given by a gamma distribution centered at age one. Rather than using a Dirac delta birth function, semelparity in this model can be achieved (or approximated) by making this distribution sufficiently narrow. In the discrete setting, the observed dichotomy between the equilibria and the boundary of the positive cone occurs due to a high co-dimension bifurcation, that is multiple eigenvalues simultaneously leave the unit circle when the extinction equilibrium destabilizes. For this continuous-time model, however, a single real eigenvalue crosses into the right-hand plane, resulting in only a branch of positive equilibria bifurcating from the extinction equilibrium. Nevertheless, we show that dynamics similar to the discrete model can occur at a nearby secondary Hopf bifurcation. In Section 2 we apply the linear chain trick to derive a continuoustime structured semelparous model. We analyze the extinction equilibrium in Section 3 and in Section 4.1 we establish the stability of the

lim ρ (t , a) = 0, a →∞

the death function is independent of age, and the birth function is separable. The assumption that the death function is independent of age is a mathematical assumption that is made in order to apply the linear chain trick to model (1) [20,21]. Biologically, this assumption may be reasonable if age-related death is negligible compared to death from hazards such as predation [28]. In addition, we assume that, as a function of age, the birth function is given by a gamma distribution centered at a = 1 whose width is determined by the parameter n ∈ N ,

n (na)ne−na. n!

(2)

This is a generalization of the nonlinearities used by Gurtin and MacCamy [19] and was given in [21] for an iteroparous species. Here larger values of n correspond to narrower distributions and, as n goes to infinity, this function converges to the Dirac delta function. Therefore, for sufficiently large n, we view this distribution as a description of a semelparous species. Since reproduction is concentrated at age one, if cycles are present in this model, we expect the cycles to have a period close to one. Next we proceed formally by defining a new set of n + 2 variables ∞

P0 (t ) =

∫0

Pi (t ) =

n (i − 1)!

ρ (t , a) da,

∫0

∞

(na)i − 1e−naρ (t , a) da for i = 1, …, n + 1.

(3)

As before, P0 denotes the total population. Meanwhile the Pi’s represent weighted measures of the population size where, as i increases, the individuals near age one (which we can interpret as the reproducing adults) are counted more heavily. Therefore, we interpret Pn + 1 as a measure (though not an exact count) of the reproducing class R and P0 − Pn + 1 as a measure of the non-reproducing class N. Given these interpretations, we assume density dependence on these two variables. The birth and death functions are now given by

β (P0, Pn + 1)

n (na)ne−na, n!

μ (P0, Pn + 1).

(4)

Note that, by the definition of the birth rate, technically individuals at any age can reproduce. However, for large n, reproduction at ages not near one is highly unlikely. Using the variables defined in (3), we reduce the McKendrick model (1) to a n + 2 -dimensional system of differential equations. This is accomplished by first integrating Eq. (1b) over the interval 0 < a < ∞ to get an equation for P0′ (t ) . Next, multiply Eq. (1b) by (na)ie−na for i = 0, …, n and, for each i, integrate the new equation over the interval 0 < a < ∞. After an application of integration by parts, an equation for Pi′ (t ) is obtained. This process results in the following system of ordinary differential equations as a model for the population dynamics,

2

dP0 = β0 (n) β (P0, Pn + 1) Pn + 1 − μ (P0, Pn + 1) P0, dt

(5a)

dP1 = nβ0 (n) β (P0, Pn + 1) Pn + 1 − (μ (P0, Pn + 1) + n) P1, dt

(5b)

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

dPi = nPi − 1 − (μ (P0, Pn + 1) + n) Pi dt

for i = 2, …, n + 1,

0 0 0 ⎛− μ 0 0 0 ⎜ 0 −(μ0 + n) ⎜ 0 ( ) 0 n μ n − + 0 J0 = ⎜ 0 0 n −(μ0 + n) ⎜ 0 0 ⎜⎜ 0 0 0 0 ⎝ 0

(5c)

where

β0 (n) =

bnn + 1 , n!

β (0, 0) = 1,

μ (0, 0) = μ0 ,

where we have used superscript zero to denote evaluation at the extinction equilibrium. Matrix (8) has characteristic polynomial

(−1)n (λ + μ0 )[(λ + n + μ0 )n + 1 − nn + 1β0 (n)]. The eigenvalues of matrix (8) are λ = −μ0 and the n + 1 roots 1 2πk ⎞ λk = −(n + μ0 ) + (β0 (n)) n + 1 n exp ⎛ i k = 0, …, n. n ⎝ +1 ⎠

Remark 1. We assume model (5) can be used to describe a semelparous species if n is sufficiently large. In addition, death soon after reproduction can be accomplished by assuming µ0 ≈ 1.

β (n ) nβ0 (n) n2β0 (n) nnβ0 (n) v (λk , β0 (n)) = col ⎛⎜ 0 , , , …, , 1⎞⎟ 2 (λ k + n + μ 0 ) n ⎠ ⎝ λ k + μ 0 λ k + n + μ 0 (λ k + n + μ 0 ) , λ k + n + μ 0 (λ k + n + μ 0 ) 2 (λ k + n + μ 0 ) n ⎞ w (λk , β0 (n)) = ⎛⎜0, , , …, , 1⎟, 2β (n ) nβ n n nnβ0 (n) ( ) 0 0 ⎝ ⎠

(10) respectively. These eigenvalues lie on a circle in the complex plane with center 1 − (n + μ0 ) and radius (β0 (n)) n + 1 n . As n increases, the center of this circle approaches negative infinity and the radius approaches positive infinity. The extinction equilibrium destabilizes when the single real root λ 0 = −(n + μ0 ) + (β0 (n))n + 1n increases through zero into the righthand plane. This occurs at the critical value β0 (n) = β0c (n) . At this critical value, the other n eigenvalues are given by

3. The extinction equilibrium We apply local bifurcation analysis to model (5) to determine the dynamics of this system. We first examine the stability of the extinction equilibrium by applying the linearization principle [17].

2πk ⎞ ⎞ λk = −(μ0 + n) ⎛1 − exp ⎛ i , k = 1, …, n, ⎝ n + 1 ⎠⎠ ⎝ ⎜

Theorem 3.1. Define

Proof. The Jacobian of (5) is

0 ⋯ 0

0 ⋯ 0 0 ⋯ 0 0 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 0

0 β0 (n) β ⎞ 0 nβ0 (n) β ⎟ ⎟ ⋯ 0 ⎟ ⋯ 0 ⎟ ⋱ 0 ⎟ n −(μ + n) ⎠

β0 (n) ∂ Pn + 1 βPn + 1 ⎞ ⎟ − ∂ Pn + 1 μP0 ⎟ nβ0 (n) ∂ Pn + 1 βPn + 1⎟ ⎟ − ∂ Pn + 1 μP1 ⎟, −∂ Pn + 1 μP2 ⎟ −∂ Pn + 1 μP3 ⎟ ⎟ ⋮ ⎟⎟ −∂ Pn + 1 μPn ⎠

(11)

In the discrete case, multiple eigenvalues simultaneously leave the unit circle when the extinction equilibrium destabilizes, resulting in both a branch of positive equilibria and a branch of synchronous cycles. This does not occur for model (5). Since a single real eigenvalue crosses into the right-half plane, only a branch of positive equilibria bifurcates from the extinction equilibrium. However at the critical value β0c (n), as n increases, the circle on which the λk eigenvalues lie gets closer to the imaginary axis. As a result, for a semelparous species with large n, a Hopf bifurcation may occur shortly after the extinction equilibrium destabilizes resulting in a branch of periodic solutions. This possibility is examined in Sections 4.2 and 5. The results in this paper generalize a result from [5], in which a secondary Hopf bifurcation was shown to occur for model (1b) with (4) and the particular nonlinearities

The extinction equilibrium of (5) is locally asymptotically stable for β0 (n) < β0c (n) and unstable for β0 (n) > β0c (n) .

⎛ β0 ∂ P0 βPn + 1 − ∂ P0 μP0 ⎜ ⎜ ⎜ nβ0 ∂ P0 βPn + 1 − ∂ P0 μP1 ⎜ +⎜ −∂ P0 μP2 ⎜ −∂ P0 μP3 ⎜ ⎜ ⋮ ⎜⎜ −∂ P0 μPn ⎝

⎟

and all have negative real parts. Therefore, by the linearization principle [17], the extinction equilibrium is stable for β0 (n) < β0c (n) and unstable for β0 (n) > β0c (n) . □

(6)

⋯ ⋯ 0 0 ⋱ 0

(9)

These eigenvalues have corresponding right and left eigenvectors

Infinite dimensional structured population models that can be reduced to a finite dimensional system of ordinary differential equations are often called linear chain trickable. These systems are particularly useful since it has been established in a general setting that the dynamics of the finite system faithfully represent the dynamics of the original model, see for instance [21,28]. Given a known solution to model (5), the density ρ(t, a) can be recovered by integrating (1b) along characteristics and applying the boundary and initial conditions defined in (1) [21].

0 0 0 ⎛− μ 0 0 ⎜ 0 −(μ + n) ⎜ 0 − + ( ) 0 n μ n J = ⎜ −(μ + n) 0 0 n ⎜ 0 0 0 ⎜ 0 0 0 ⎝ 0

0 β0 (n) ⎞ 0 nβ0 (n) ⎟ ⎟ 0 ⋯ ⎟, 0 ⋯ ⎟ 0 ⋱ ⎟ n −(μ0 + n) ⎟⎠ (8)

and where we have dropped the time dependence for notational convenience. By definition, 1/µ0 is the average lifespan in the absence of density effects. Since reproduction is assumed to occur only on a small interval around age one, the assumption µ0 > > 1 means that, on average, most individuals do not survive to become reproductive. If this is the case, it is biologically unlikely for cycles to occur. For many semelparous species, reproduction is quickly followed by death. This can be accomplished by assuming µ0 near one.

μ n+1 β0c (n) ≜ ⎛1 + 0 ⎞ . n⎠ ⎝

⋯ ⋯ 0 0 ⋱ 0

β (P0, Pn + 1) = max{1 − Pn + 1, 0},

μ (P0, Pn + 1) = μ0 ,

when n > 1. 4. The positive equilibria 4.1. The stability of positive equilibria As β0(n) increases past the critical value β0c (n), a branch of positive equilibria bifurcates from the extinction equilibrium [6]. Setting β0 (n) = β0c (n)(1 + ɛ), we consider the equilibria as parametrized by ε and parameterize the bifurcating branch by (β0(n)(ε), Pi(ε)) where (β0 (n)(0), Pi (0)) = (β0c (n), 0) . This method requires the following

(7)

where we use the shorthand μ = μ (P0, Pn + 1) and β = β (P0, Pn + 1) . When evaluated at the extinction equilibrium, the Jacobian becomes 3

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

where λ 0 (0) = 0 and λk(0), given in (11), is obtained by setting β0 (n) = β0c (n) in (9). Standard Lyapunov-Schmidt expansion gives the calculation

smoothness assumptions, A1. For the open set Ω ⊂ R with R+ ⊂ Ω, assume

β (P0, Pn + 1), μ (P0, Pn + 1) ∈ C 2 (Ω × Ω → R+1 ), β (P0, Pn + 1) ≤ βmax < ∞, 0 ≤ μ (P0, Pn + 1).

λk′ (0) =

In addition, the stability of the positive equilibria relies of the term

a+ (n) ≜ (n + μ0 ) [

μ0 ∂0Pn + 1 β

+

+

β0c (n) ∂0P0 β

] − (n + 1) [

μ0 ∂0Pn + 1 μ

].

β0c (n) ∂0P0 μ

wT (λk (0), β0c (n)) J ′ (0) v (λk (0), β0c (n))

Theorem 4.1. Assume A1. At β0 (n) = a branch of positive equilibria bifurcates from the extinction equilibrium of (5) if a+ (n) ≠ 0 .

B (0) =

Pi (ɛ) = 0 + Pi′ (0)ɛ + O (ɛ2).

⎛ β0 (n) Pn′+1 ∂ P0 β ⎜ ⎜nβ0c (n) Pn′+1 ∂0P0 β ⎜ ⎜ 0 ⋮ ⎜ 0 ⎝ ′∂0

The equilibria equations for model (5) are

(13b)

0 = nPi − 1 − (μ (P0, Pn + 1) + n) Pi for i = 2, …, n + 1.

(13c)

M (0) =

⎛ −P0 0P0 μ ⎜ −P1′∂ P0 μ 0 ⎜ −P2′∂ P0 μ ⋮ ⎜ −Pn′+ 1 ∂0P0 μ ⎝

β0c (n) + 2β0c (n) Pn′+ 1 ∂0Pn + 1 β

0 ⋯ 0

⎞ ⎟ nβ0c (n) + 2nβ0c (n) Pn′+ 1 ∂0Pn + 1 β ⎟ , ⎟ + nβ0c (n) P0′∂0P0 β ⎟ 0 ⋮ ⎟ 0 ⎠ + β0c (n) P0′∂0P0 β

0 ⋯ 0 0 ⋯ 0 ⋮ ⋮ ⋮ 0 ⋯ 0

0 ⋯ 0 0 ⋯ 0 0 ⋯ 0 ⋮ ⋮ ⋮ 0 ⋯ 0

−P0′∂0Pn + 1 μ

⎞ ⎟ −P2′∂0Pn + 1 μ , ⎟ ⋮ ⎟ −Pn′+ 1 ∂0Pn + 1 μ ⎠ −P1′∂0Pn + 1 μ

and D (0) = diag (−Pn′+ 1 ∂0Pn + 1 μ − P0′∂0P0 μ) is a diagonal matrix whose diagonal components are given by the term − Pn′+ 1 ∂0Pn + 1 μ − P0′∂0P0 μ . Note that for any eigenvalue the following simplifications hold

Solving Eq. (13b) for P1 and then Eq. (13c) for Pi, we get i

n ⎞⎟ Pn + 1, 1 ≤ i ≤ n. Pi = β0 (n) β (P0, Pn + 1) ⎛⎜ ⎝ n + μ (P0, Pn + 1) ⎠

wT (λk (0), β0c (n)) v (λk (0), β0c (n)) = n + 1,

We substitute the solution for Pn into Eq. (13c) for i = n + 1. We can now reduce the system to two equilibria equations in terms of P0 and Pn + 1. Setting β0 (n) = β0c (n)(1 + ɛ) and parameterizing the equilibria by ε, we are left with solving the following equations

0 =

0

c

Proof. To determine the existence and direction of bifurcation of positive equilibria, we calculate the Taylor expansion of (β0(n)(ε), Pi(ε)) to order one in ε at ɛ = 0,

0 = nβ0 (n) β (P0, Pn + 1) Pn + 1 − (μ (P0, Pn + 1) + n) P1,

(17)

where

(a) If a+ (n) > 0, then the branch of positive equilibria bifurcates backward and is unstable. (b) If a+ (n) < 0, then the branch of positive equilibria bifurcates forward and is stable for β0 (n)⪆β0c (n) .

(13a)

(16)

J ′ (0) = B (0) + M (0) + D (0),

β0c (n),

0 = β0 (n) β (P0, Pn + 1) Pn + 1 − μ (P0, Pn + 1) P0,

,

where J′(0) is the ε-derivative of the Jacobian (7) with β0 (n) = β0c (n)(1 + ɛ) evaluated at ɛ = 0 . See [6] for details on this calculation. The matrix J′(0) can be decomposed into three matrices,

(12)

This term is a weighted sum of the low-density sensitivities representing a measure of the total competition in the population.

β0c (n)(1

wT (λk (0), β0c (n)) v (λk (0), β0c (n))

wT (λk (0), β0c (n)) D (0) v (λk (0), β0c (n)) n+1

(18a)

= − (Pn′ + 1 (0) ∂0Pn + 1 μ + P0′ (0) ∂0P0 μ) =

+ ɛ) β (P0 (ɛ), Pn + 1 (ɛ)) Pn + 1 (ɛ) − μ (P0 (ɛ), Pn + 1 (ɛ)) P0 (ɛ),

(n + μ0 ) (μ0 ∂0Pn+1 μ + β0c (n) ∂0P0 μ), a+ (n)

(18b)

(14a) n

n ⎞ − (μ (P0 (ɛ), Pn + 1 (ɛ)) + n). 0 = β0c (n)(1 + ɛ) ⎜⎛ ⎟ ⎝ n + μ (P0 (ɛ), Pn + 1 (ɛ)) ⎠ (14b)

wT (λk (0), β0c (n)) M (0) v (λk (0), β0c (n)) n+1

n+1

(λk (0),

∑ P′j (0) wj +1 j=1

β0c (n)),

(18c)

P0 (ɛ) = − β0c (n)

where wj + 1 denotes the j + 1 component of vector w given by Eq. (10). Since the destabilization of the extinction equilibrium occurs when a single real eigenvalue λ0 increases past zero, the real parts of all other eigenvalues are negative and remain negative for ε ≳ 0. Therefore, the stability of the positive equilibria depends only on the sign of e(λ 0′ (0)) . We calculate λ 0′ (0) using (17) and (10) in (16). From the decomposition of J′(0) we obtain

(15)

Therefore, positive equilibria bifurcate backwards (ε < 0) if a+ (n) > 0 and forwards (ε > 0) if a+ (n) < 0 . To determine the local stability of the bifurcating equilibria, we consider the eigenvalues of the Jacobian (7) evaluated at the positive equilibrium (15) for ε ≈ 0. These eigenvalues are given by the Taylor expansion

λk (ɛ) = λk (0) + λk′ (0)ɛ +

c −1 ⎛ β0 (n) 0 ∂P μ ⎜ n + 1 ⎝ λk (0) μ0 0

+ ∂0Pn + 1 μ ⎟⎞ ⎠

We differentiate the equilibria equations (14) with respect to ε, evaluate the result at ɛ = 0, and solve for Pi′ (0) . This results in the equilibria expansions

n + μ0 ɛ + O (ɛ2), a+ (n) n + μ0 ɛ + O (ɛ2), Pn + 1 (ɛ) = − μ0 a+ (n) n + μ0 ni ɛ + O (ɛ2) for 1 ≤ i ≤ n. Pi (ɛ) = − β0c (n) μ0 a+ (n) (n + μ0 )i

=

λ 0′ (0) = −

n + μ0 ɛ, n+1

where we have used (18) to obtain

O (ɛ2), 4

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

wT (λk (0), β0c (n)) B (0) v (λk (0), β0c (n))

λ1′ (0) =

n+1 −(n + μ0 ) [(n + μ0 ) (μ0 ∂0Pn + 1 β + β0c (n) ∂0P0 β )+(n = (n + 1) a+ (n) + 1) (μ0 ∂0Pn + 1 μ + β0c (n) ∂0P0 μ) ]

∑ P′j (0) wj+1 (λ 0 (0), β0c (n)) = − j=1

+ ∂0Pn + 1 μ ⎟⎞ ⎠

∑ P′j (0) wj+1 (λ1 (0), β0c (n))

= − μ0

j=1

(n + μ0 ) ⎡ 1 a+ (n) ⎢ ⎣ n

j

λ1 (0) + n + μ0 ⎞ ⎤ ⎟ ⎥ n + μ0 j=1 ⎝ ⎠⎦ (n + μ0 ) n 2πj = − μ0 ∑ exp ⎛ n + 1 i ⎞. a+ (n) j = 0 ⎝ ⎠ +

∑⎛

⎜

Since this term is a multiple of the sum of the n + 1 roots of unity, it equals zero. Therefore, the third term disappears and we obtain

For large n, the first complex pair of eigenvalues associated with the extinction equilibrium (corresponding to k = 1 and k = n ) are close to the imaginary axis. This suggests, for sufficiently large n, that as β0(n) increases, a pair of complex eigenvalues associated with the positive equilibria might cross into the right-hand plane resulting in a Hopf bifurcation. Ideally, we would calculate the eigenvalues of the Jacobian (7) evaluated at the positive equilibrium (15). However, there is no closed form for the zeros of the characteristic equation of the Jacobian when evaluated at a positive equilibrium. Instead, we calculate the order ε approximate of the eigenvalues at the extinction equilibrium,

wT (λ1 (0), β0c (n)) B (0) v (λ1 (0), β0c (n)) ⎞  m(λ1′ (0)) =  m ⎜⎛ ⎟ n+1 ⎠ ⎝ =

s (n) ⎡−(μ (n + μ ) ∂0 β + (n + 1) (μ ∂0 μ 0 0 Pn + 1 0 Pn + 1 a+ (n)(n + 1) ⎢ ⎣ μ0 β0c (n) n (n + μ0 ) ∂0P0 β ⎤ + β0c (n) ∂0P0 μ)) + , (μ0 + r (n))2 + s (n)2) ⎥ ⎦ (21) T ⎛ w (λ1 (0),

e(λ1′ (0)) = e ⎜ ⎝

β0 (n) = β0c (n)(1 + ɛ).

β0c (n)) B (0) v (λ1 (0),

β0c (n)) ⎞ ⎟

n+1

⎠

(n + μ0 ) + (μ0 ∂0Pn+1 μ + β0c (n) ∂0P0 μ), a+ (n)

This estimate of λi(ε) is appropriate for sufficiently small ε. In particular, the first eigenvalues to cross into the right-hand plane will be λ1(ε) and λn(ε), so we restrict out attention to λ1(ε). Recall that e (λ1(0)) < 0 and ε > 0 for forward bifurcating positive equilibria. Therefore, if e(λ1′ (0)) > 0, then λ1(ε) has zero real part for

(22)

where the first term of (22) calculates to

−1 [(n + r (n) + μ0 ) (μ0 (n + μ0 ) ∂0Pn + 1 β (n + 1) a+ (n) + (n + 1) (μ0 ∂0Pn + 1 μ + β0c (n) ∂0P0 μ))

(19)

(μ0 + r (n)) n c 0 ⎤ ⎞ + ⎜⎛1 + ⎟ (n + μ 0 ) β (n ) μ 0 ∂ P β . 0 ⎥ 0 2 2 ( μ 0 + r (n )) + s (n ) ⎠ ⎝ ⎦

If εh is sufficiently small, then e(λ1 (ɛ)) = 0 for β0(n) not too far from β0c (n) . In particular, when e(λ1′ (0)) > 0, λ1 (ɛ) has positive real part for ε ≳ εh. Since ɛ = ɛh corresponds to the existence of a simple pair of imaginary eigenvalues at the critical value

Combining the terms of e(λ1′ (0)) and using the definitions of r(n) and s(n) given by Eq. (20), Eq. (22) can be rewritten as

e(λ1′ (0)) =

β0h (n) ≜ β0c (n)(1 + ɛh),

−(n + μ0 ) ⎡ ⎛ 2π ⎞ ⎞ (n + 1) (μ0 ∂0Pn + 1 μ −1 + cos ⎛ (n + 1) a+ (n) ⎢ + 1 ⎠⎠ n ⎝ ⎣⎝ ⎜

⎟

+ β0c (n) ∂0P0 μ)

this indicates a Hopf bifurcation. We next calculate e(λ1′ (0)) and  m(λ1′ (0)) . For ɛ = 0, λ1 (0), given in Eq. (11), has real and imaginary parts

2π ⎞ r (n) ≜ e(λ1 (0)) = −(n + μ0 ) + (n + μ0 )cos ⎛ , ⎝n + 1⎠

j=1

n+1

4.2. Destabilization of the positive equilibria

e(λ1 (0)) . e(λ1′ (0))

∑ P′j (0) wj+1 (λ1 (0), β0c (n)),

where the second term, which follows from (18), is real and negative by A2. Simplification of the sum in the third term shows

Notice that the bifurcation of positive equilibria is forward if we ∂0P0 β ≤ assume only negative density effects, 0 0 0 0, ∂ Pn + 1 β ≤ 0, ∂ P0 μ ≥ 0, ∂ Pn + 1 μ ≥ 0, and at least one partial derivative is non-zero. Since this is the case we are interested in, we make that assumption for the remainder of this paper. A2. Assume a+ (n) < 0 .

ɛh ≜ −

c 1 ⎛ β0 (n) ∂0P μ ⎜ n + 1 ⎝ λ1 (0) + μ0 0 n+1

μ0 (n + 1)(n + μ0 ) . a+ (n)

The backward bifurcating equilibria (ε < 0) are unstable since λ 0′ (0) > 0 . Meanwhile, the forward bifurcating equilibria (ε > 0) are stable since λ 0′ (0) < 0 . □

λi (ɛ) = λi (0) + λi′ (0)ɛ,

n+1 (n + μ0 ) (μ0 ∂0Pn+1 μ + β0c (n) ∂0P0 μ) + a+ (n) −

and n+1

wT (λ1 (0), β0c (n)) B (0) v (λ1 (0), β0c (n))

2π ⎞ 0 + μ0 (n + μ0 )cos ⎛ ∂ Pn + 1 β ⎝n + 1⎠ (μ0 + r (n)) n 0 ⎤ ⎞ c + ⎜⎛1 + ⎟ β (n ) μ 0 ∂ P β . 0 ⎥ 0 2 2 ( μ 0 + r (n )) + s (n ) ⎠ ⎝ ⎦

(20a)

(23)

2π ⎞ s (n) ≜  m(λ1 (0)) = (n + μ0 )sin ⎛ . ⎝n + 1⎠

Since a+ (n) < 0, the expression in the square brackets in (23) must be positive to have e(λ1′ (0)) > 0 . If e(λ1′ (0)) > 0 and λ (ɛh) = 0 for sufficiently small εh, then the positive equilibria are unstable for ε > εh. In the following sections, we use Eqs. (21) and (23) to study and interpret the nature of the bifurcation that occurs at β0h (n) . In order to relate our results to the discrete semelparous model in which the

(20b)

Using (10) and (17) in (16) we obtain

5

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

• If d < 0, then the positive equilibria are stable for ε > ε

stability conditions for the positive equilibria and synchronous 2-cycles depend on the relative levels of between-class and within-class competition, we first introduce a change of variables. To allow for the comparison of these two types of competition, we define

R ≜ Pn + 1,

• •

N ≜ P0 − Pn + 1,

to represent measures of the reproducing class R and the non-reproducing class N. Accordingly, we assume that the density terms have the following density dependence

h and unstable for ε < εh. If d > 0, then the positive equilibria are stable for ε < εh and unstable for ε > εh. If d · l1(εh) > 0, then the bifurcation of periodic solutions is backward and if d · l1(εh) < 0, then the bifurcation is forward. The periodic solutions have the opposite stability conditions as the positive equilibria, provided they exist.

In this section we consider the first two conditions for the existence of a Hopf bifurcation. Calculation of the first Lyapunov exponent in (iii) is not tractable for this model. However, we use numerical simulations in the next section to verify the existence of stable periodic solutions.

β (P0 − Pn + 1, Pn + 1) = β (N , R), μ (P0 − Pn + 1, Pn + 1) = μ (N , R). From this change of variable we have

∂0P0 β = ∂0N β , ∂0Pn + 1 β = ∂0R β − ∂0N β ,

5.1. The non-hyperbolicity condition

∂0P0 μ = ∂0N μ, ∂0Pn + 1 μ = ∂0R μ − ∂0N μ,

We first show that if e(λ1′ (0)) > 0 and n is sufficiently large, then there exists an εh ≳ 0 satisfying condition (i).

which allows us to rewrite Eq. (23) as

Proposition 5.1. Assume e(λ1′ (0)) > 0 . If n is sufficiently large, then condition (i) is satisfied. In addition, if a Hopf bifurcation occurs, the bifurcation occurs sooner (that is, for β0h (n) closer to β0c (n) ) for larger values of n.

−(n + μ0 ) ⎡ 2π ⎞ ⎞ c (n + 1) ⎛−1 + cos ⎛ [(β0 (n) e(λ1′ (0)) = (n + 1) a+ (n) ⎢ ⎝ n + 1 ⎠⎠ ⎝ ⎣ ⎜

⎟

− μ0 ) ∂0N μ + μ0 ∂0R μ]

Proof. To determine when condition (i) is satisfied we examine the leading term behavior of e(λ1′ (0)) and  m(λ1′ (0)) . Consider e(λ1′ (0)) given in Eq. (24). Since a+ (n) has order O(n), the factor outside of the square brackets in (24) has order O(1/n). Using the Taylor series for large n to approximate

2π ⎞ 0 [∂ R β − ∂0N β ] + μ0 (n + μ0 )cos ⎛ n ⎝ + 1⎠ (μ0 + r (n)) n 0 ⎤ ⎞ c + ⎜⎛1 + ⎟ β (n ) μ 0 ∂ J β , ⎥ (μ0 + r (n))2 + s (n)2 ⎠ 0 ⎝ ⎦ (24)

2π ⎞ 2π ⎞ 2π cos ⎛ ≈ 1 + O (1/ n2), sin ⎛ ≈ + O (1/ n3), n+1 ⎝n + 1⎠ ⎝n + 1⎠

where

we find that the three terms inside the square brackets of Eq. (24) have orders O(1/n), O(n) and O(n), respectively. Therefore, e(λ1′ (0)) has order O(1). In addition, for sufficiently large n, the death density dependence has vanishingly small effect on e(λ1′ (0)) . Therefore, the sign of e(λ1′ (0)) is determined by the density effects on the birth rate. In particular, for sufficiently large n, the sign of e(λ1′ (0)) is equal to the sign of

a+ (n) = (n + μ0 )[μ0 ∂0R β + (β0c (n) − μ0 ) ∂0N β ] − (n + 1)[μ0 ∂0R μ + (β0c (n) − μ0 ) ∂0N μ], and the expanded form of the factor in front of the last term in (24) is given by

1+

(μ0 + r (n)) n (μ0 + r (n))2 + s (n)2

(−n + (n + μ )cos ( ) ) n . =1+ (−n + (n + μ )cos ( ) ) + ( (n + μ )sin ( ) ) 0

0

2π n+1

2π ⎞ 0 [∂ R β − ∂0N β ] μ0 (n + μ0 )cos ⎛ ⎝n + 1⎠ (μ0 + r (n)) n 0 ⎞ c + ⎜⎛1 + ⎟ β (n ) μ 0 ∂ N β . 0 2 + s (n ) 2 ( ( )) μ r n + 0 ⎠ ⎝

2π n+1

2

0

2π n+1

(25)

2

(26)

Recall that, using the order ε Taylor expansion of e(λ1(ε)), we may solve e(λ1 (ɛ)) = 0 to get εh as defined in Eq. (19). Since e (λ1(0)) has order O(1/n) and e(λ1′ (0)) has order O(1), as n approaches infinity, εh approaches zero at a rate of 1/n. Additionally, since e(λ1(0)) < 0, εh is positive if e(λ1′ (0)) > 0 . Therefore, if e(λ1′ (0)) > 0 and n is sufficiently large, then there exists εh ≳ 0 such that u (ɛh) = 0 . To satisfy condition (i), it remains to verify that v(ɛh) ≠ 0, that is s (n) +  m(λ1′ (0))ɛh ≠ 0 . Consider  m(λ1′ (0)) given in Eq. (21). The term outside of the square brackets in (21) has order O(1/n2) and the terms inside the square brackets have orders O(1) and O(n2), respectively. Therefore  m(λ1′ (0)) has order O(1). Since s(n) has order O(1) and εh has order O(1/n), v(ɛh) has order O(1) and approaches 2π as n goes to infinity. Hence, v(ɛh) ≠ 0 for sufficiently large n. □

5. Hopf bifurcations for semelparous species Recall that, in our modeling methodology, large n corresponds to a semelparous species. In this section, we restrict our focus to the case of large n and examine the conditions under which a Hopf bifurcation might occur when the positive equilibria destabilize. In order to establish the existence and stability of a Hopf bifurcation, there must exist an εh ≳ 0 such that the following conditions are satisfied [30]:

λ1, n (ɛ) = u (ɛ) ± (i) Non-hyperbolicity condition: Given iv (ɛ), u (ɛh) = 0, v (ɛh) ≠ 0 and all other eigenvalues have negative real parts. du (ɛ) (ii) Transversality condition: dɛ ɛ = ɛh = d ≠ 0 . (iii) Genericity condition: The first Lyapunov exponent is non-zero, l1(εh) ≠ 0.

Remark 2. Assume a Hopf bifurcation occurs for εh ≳ 0 resulting in 2π stable periodic cycles. Then v (ɛh) ∼ (n + μ0 ) n + 1 and the cycles have period

Conditions (i) and (ii) establish the existence of a unique branch of periodic solutions that bifurcates from the positive equilibria with period 2π/v(ɛh). Meanwhile, condition (iii) is needed to determine the direction of bifurcation and stability of the periodic solutions. If these three conditions are satisfied, then a Hopf bifurcation occurs at εh that has the following properties [30]:

n+1 n + μ0

which converges to one as n goes to infinity.

5.2. The transversality condition Next we calculate the sign of d = e(λ1′ (0)) for large n to verify conditions (i) and (ii) are satisfied. 6

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

Proposition 5.2. For sufficiently large n, e(λ1′ (0)) is non-zero and satisfies condition (ii). In particular, for sufficiently large n and small µ0, if ∂0R β is sufficiently large relative to ∂0N β , then e(λ1′ (0)) > 0 and the positive equilibria are unstable for ε ≳ εh. If, in addition l1(0) < 0, then periodic cycles are stable for ε ≳ εh.

Eq. (10). At εh,

2π ⎞ λ1 (ɛh) = (n + μ0 )sin ⎛ +  m(λ1′ (0))ɛh ⎝n + 1⎠ which, for large n, is dominated by λ1(0) and, using (25), can be approximated as

Proof. Recall that we are assuming only negative density effects which means that a+ (n) < 0 and the factor outside of the square brackets in (24) is positive. In addition, since β0c (n) converges to e μ0 from above, β0c (n) − μ0 is positive. Therefore, the coefficient in front of the first term of (24) is negative and the coefficients in front of the second and third terms are positive. Overall, the first and third terms contribute negatively to e(λ1′ (0)) . Meanwhile, the sign of the second term depends on the sign of ∂0R β − ∂0N β . It follows that a necessary condition for e(λ1′ (0)) > 0 is ∂0R β − ∂0N β > 0. Next we verify that it is possible for (26) to be positive. Using the approximations in (25) and β0c (n) ∼ e μ0, we get

1+

λ1 (ɛh) ∼ 2πi. Normalizing the first entry of v(λ1 (ɛh), β0c (n)), we have

v0 (λ1 (ɛh), β0c (n)) = 1 and vn + 1 (λ1 (ɛh), β0c (n)) = Using the approximation β0c (n) ∼ e μ0, we obtain

P0 (t ) = P0 (ɛ) + δ cos(2πt ), 1 Pn + 1 (t ) = Pn + 1 (ɛ) + δ μ [μ0 cos(2πt ) − 2π sin(2πt )]. e 0

(μ0 + r (n)) n μ0 n (n + 1)2 , ∼1+ 2 2 2 (μ0 + r (n)) + s (n) μ0 (n + 1)2 + 4π 2 (n + μ0 )2

Finally, we rewrite these in terms of the variables R and N using the trigonometric identity

and can approximate (26) by

⎞ 0 μ0 (n + μ0 ) ∂0R β + ⎛⎜−μ0 (n + μ0 ) + 1 + 2 ⎟ ∂ N β. μ0 (n + 1)2 + 4π 2 (n + μ0 )2 ⎠ ⎝

to get

For sufficiently large n, this term grows like

R (t ) = R (ɛ) + δ

e μ0

μ ⎡ ⎤ − 1⎞⎟ ∂0N β ⎥. nμ0 ⎢∂0R β + ⎜⎛ 2 0 2 + μ 4 π 0 ⎝ ⎠ ⎣ ⎦

μ02 + 4π 2

N (t ) = N (ɛ) + δ

Therefore, for sufficiently large n, a necessary condition for e(λ1′ (0)) > 0 is

μ02 + 4π 2

b a2 + b2 cos ⎡2πt − tan−1 ⎛ ⎞ ⎤ ⎢ a ⎝ ⎠⎥ ⎣ ⎦

a cos(2πt ) + b sin(2πt ) =

μ02 e μ0 n (n + 1)n

μ0 e μ0

μ0 + 2πi . β0c (n)

−1 ⎛ 2π ⎞ ⎤ cos ⎡ ⎢2πt + tan ⎜ μ ⎟ ⎥, ⎝ 0 ⎠⎦ ⎣ − μ0 )2 + 4π 2 2π ⎞ ⎤ ⎡ cos ⎢2πt − tan−1 ⎜⎛ μ ⎟ . ⎥ 0 − μ e μ0 e 0 ⎠⎦ ⎝ ⎣

e μ0 (e μ0

These equations can be used to examine whether or not the reproducing and non-reproducing cycles are out of sync.

− 1 < 0. (27)

Remark 4. The difference in phase shifts defined by

□

Δϕ (μ0 ) ≜ ϕN − ϕR =

We note that Eq. (27) holds for µ0 ≲ 2.82. Numerically, we can verify that n does not need to be very large to get e(λ1′ (0)) > 0 . For instance, choosing μ0 = 0.8, ∂0N β = −4, ∂0N β = −1, and ∂0N μ = ∂0R μ = 0, we have e(λ1′ (0)) > 0 for n > 4.

1 ⎡ −1 ⎛ 2π ⎞ 2π ⎞ ⎤ tan ⎜ ⎟ + tan−1 ⎜⎛ μ ⎟ ⎥ 0 − μ 2π ⎢ μ e 0 0 ⎠⎦ ⎝ ⎠ ⎝ ⎣

is a monotonically decreasing function in µ0 that gives a measure of how out of sync are the periodic cycles R (t ) and N (t ). Since these cycles have period one, Δϕ (μ0 ) = 0.5 corresponds to being completely out-of-phase.

Remark 3. A necessary condition for e(λ1′ (0)) > 0, assuming only negative feedback and large n, is

Fig. 1(a) shows Δϕ(µ0) for 0 ≤ µ0 ≤ 3. For μ0 = 1, this difference is equal to 0.43 radians. Fig. 1(b) shows the ratio of the amplitudes of N (t ) and R (t ), AN / AR , for 0 ≤ µ0 ≤ 3. For µ0 near one, this ratio is close to one which means that the cycles have approximately the same amplitude.

∂0R β − ∂0N β > 0. This condition means that the effect of non-reproducers on adult reproduction (between-class competition) is greater than the effect of reproducers on adult reproduction (within-class competition).

6. A numerical study 5.3. Calculating periodic orbits near the Hopf bifurcation

where δ ≈ 0 and

In Section 5 we showed that, if 0 < µ0 < 2.82 and n is sufficiently large, then it is possible for the stable positive equilibria to destabilize via a Hopf bifurcation for β0(n) near β0c (n) . Whether or not a Hopf bifurcation occurs depends on the density effects on the birth rate. In particular, the between-class competition ∂0N β must be greater than the within-class competition ∂0R β . Additionally, when a Hopf bifurcation occurs, the length of the interval over which the positive equilibria are stable decreases as n increases. In this section we demonstrate these results and study numerically whether or not the resulting periodic cycles are stable. We also verify that the non-reproducing and reproducing classes are out of sync in these periodic cycles, as occurs in the discrete setting. In what follows we use density terms

hi (t ) ≜ e(vi (λ1 (ɛh), β0c (n)))cos(2πt ) −  m(vi (λ1 (ɛh), β0c (n)))sin(2πt ).

β (N , R) = exp(−c1 N − c2 R),

Here we use vi to denote the ith entry of the eigenvector v given in

For large n, the density effects on the death rate have minimal effect on

Recall that, given sufficiently strong between-class competition, the discrete Leslie model can exhibit stable synchronous cycles in which the reproducing and non-reproducing classes are separated temporally. The equivalent to a synchronous cycle in the continuous model (5) is a stable periodic cycle in which the reproducing and non-reproducing classes are out of phase. Here we verify that this occurs for large n and µ0 near one. Periodic orbits near a Hopf bifurcation can be approximated by Seydel [30]

Pi (t ) = Pi (ɛ) + δhi (t ),

7

μ (N , R) = exp(c3 N + c4 R),

ci ≥ 0.

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

Fig. 1. (a) Shown is the difference in phase shifts Δϕ (µ0) for 0 ≤ µ0 ≤ 3. The periodic cycles for R and N are completely out of sync if Δϕ (μ 0 ) = 0.5. (b) Shown is the ratio of the amplitudes of the cycles AN/AR. These cycles have the same amplitude when AN / AR = 1.

Fig. 2. Shown are the time series of R and N when n = 50 . Plots (c) and (d) show the transients for 10 initial time steps and the final attractor after 300 or 50 time steps. For all graphs, μ 0 = 0.8, c3 = 0, c4 = 0, and

β0c (n) ≈ 2.24 . Additional parameter values are row 1: c1 = −1 and c2 = −4, row 2: c1 = − 4 c2 = −1, and column 1: and column 2: β0(50) ≈ 11.23, β0(50) ≈ 22.47. In (a) and (b), a positive equilibrium is stable for both values of β0(50). Meanwhile, in (c) an equilibrium is stable for small values of β0(50) but destabilizes in (d) to a periodic cycle for larger values of β0(50).

the dynamics. Therefore, we first consider the case where density effects only occur in the birth rate (c3 = c4 = 0 ). Choosing μ0 = 0.8 and n = 50, Fig. 2 shows the time series for R and N and demonstrates that strong between-class competition can cause a Hopf bifurcation resulting in stable periodic cycles. In this figure, if ∂0R β is sufficiently larger than ∂0N β , then a positive equilibrium is stable in a neighborhood of the extinction bifurcation. This is shown in Fig. 2(a) and (b) for two values of β0(50). However, if ∂0R β is sufficiently smaller than ∂0N β , then the equilibrium is stable near the extinction bifurcation, shown in Fig. 2(c), but destabilizes into a cycle further away, shown in Fig. 2(d). Note that the period of these cycles is close to one (approximately 0.94) and that R and N are out of sync. Here the time between the maximum of R and the maximum of N is approximately 0.41 (completely out of sync would be the half of the period or 0.47). Biologically, we would expect many semelparous species to have µ0 near one since adults typically die shortly after reproducing. In Fig. 3, we show that µ0 needs to be near 1 in order for cycles to occur. Increasing µ0 close to the endpoints of the interval 0 < µ0 < 2.82 can remove the Hopf bifurcation that occurs in Fig. 2. In Fig. 3, we consider the case of μ0 = 1.8 (row 1) and μ0 = 0.2 (row 2). Choosing the same β0(50) value

as in Fig. 2(d), we see that a Hopf bifurcation no longer occurs (column 1). However, increasing β0(50) to a value further from β0c (n), results in a stable cycle for μ0 = 0.2 but not μ0 = 1.8 (column 2). Next, we consider the effect of including a density dependent death rate. In Fig. 4 we use c3 = c4 = 0.5 to show that density dependent death rates do not promote population cycles. Keeping n = 50 and the other parameter values as in Fig. 2(d), a Hopf bifurcation no longer occurs, as shown in Fig. 4(a). This is due to the negative contribution the death density terms make to e(λ1′ (0)) . However, if we increase n to n = 100, the density effects from the death rate have less influence on the dynamics, resulting in a stable cycle, as shown in Fig. 4(b). Finally, we show that the distance between the Hopf bifurcation and the primary bifurcation, as measured by the difference β0h (n) − β0c (n), decreases with increasing n. Fig. 5 shows the bifurcation diagrams for n = 15, n = 50, and n = 100 . Choosing the same parameter values as Fig. 2(d), we see in Fig. 5(a) that for n = 15 no Hopf bifurcation occurs. However, if we increase n, not only does a Hopf bifurcation occur, but it occurs sooner (for smaller values of β0(n)) for larger n, as shown in Fig. 5(b) and (c).

8

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

Fig. 3. Shown are the time series of R and N when n = 50 . For all graphs, c1 = −4, c2 = −1, c3 = 0, c4 = 0 . Plot (d) shows the transients for 40 initial time steps and the final attractor after 150 time steps. Additional parameter values are row 1: μ 0 = 1.8 (resulting in β0c (n) ≈ 6.07 ), row 2: μ 0 = 0.2 (resulting in β0c (n) ≈ 1.22 ), column 1: β0(50) ≈ 11.23, and column 2: β0(50) ≈ 67.40. When µ0 is closer to the endpoints of the interval 0 < µ0 < 2.82, a greater difference in the competition terms c1 and c2 is required for a Hopf bifurcation to occur. In (a), (b) and (c) a Hopf bifurcation does not occur. However, in (d), a Hopf bifurcation occurs for a value of β0(50) further away from β0c (n) .

7. Concluding remarks

(i) For large n and small µ0, 0 < µ0 < 2.82, the positive equilibria destabilize via a Hopf bifurcation if the amount of between class competition is sufficiently greater than the amount of within class competition. (ii) Numerical simulations show that the Hopf bifurcation can result in stable periodic cycles. Cycles are easier to obtain for µ0 closer to one. In addition, the two classes are out of sync and the cycles have period near one (the maturation period). (iii) The destabilization of the positive equilibria occurs sooner (for smaller values of β0(n)), for larger n, that is when the species exhibit greater semelparity as defined by the width of the birth function. (iv) For large n, the effect of density on the death rate has vanishingly small effect on system dynamics near the extinction bifurcation.

In this paper we present model (5) as one possible continuous-time model for a semelparous population. However, other model assumptions could be used to derive a continuous-time semelparous model. For instance, one may wish to use a different birth function such as a Dirac delta function to represent semelparity. We explored this possibility initially but it proved to be intractable. The birth function used in this paper appears to be an appropriate substitute given that it converges to a Dirac function as n goes to infinity. In addition, there are other possible definitions for the reproducing class. For example, for large n, Pn is not very different from Pn + 1. Therefore, one may wish to represent the reproducing class as a weighted sum of the Pk through Pn + 1 classes. This also proves difficult to analyze analytically but numerical simulations (not shown) using a weighted sum of the last two classes suggest that this choice of density dependence gives qualitatively similar results. The main results for model (5) are summarized by the following:

The combination of (i)–(iii) establishes, in a continuous model, dynamics similar to the synchronous 2-cycles that occur in the discretetime, two-stage semelparous Leslie model. The last result (iv) does not Fig. 4. Shown are the time series of R and N when there are non-zero death density effects. Both plots show the transients for 10 initial time steps and the final attractor after 150 or 50 time steps. For both graphs, c1 = −4, c2 = −1, c3 = 0.5, c4 = 0.5, μ 0 = 0.8. Additional parameter values are (a) n = 50 (resulting in β0c (n) ≈ 2.24 ) and β0 (50) = 22.47, and (b) n = 100 (resulting in β0c (n) ≈ 2.24 ) and

β0 (50) = 22.36 . For small n, the death density effects result in a stable equilibrium, as shown in (a). For sufficiently large n, the effects of the death density terms becomes negligible, resulting in a stable periodic cycle, as shown in (b).

9

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

Fig. 5. Shown are the bifurcation diagrams for n = 15, n = 50, and n = 100 . Additional parameter values are μ 0 = 0.8, c1 = −4, c2 = −1, c3 = 0, c4 = 0. In (a), n = 15 results in stable equilibria for the range of β0(n) values considered. In (b) and (c), the positive equilibria destabilize via a Hopf bifurcation for sufficiently large n. Additionally, this bifurcation occurs sooner (that is the interval of β0(n) values over which the positive equilibria are stable is shorter) for larger n. This can be seen by comparing (b) and (c).

occur for the discrete semelparous model. However, note that in the continuous model, every age has the same death rate. Therefore µ represents both a between-class and within-class term and does not have a direct equivalent in discrete models. Overall, model (5) provides a continuous-time counterpart to the two-stage, discrete-time semelparous model. Though it does not exhibit equivalent dynamics to the discrete-time models, it does share many features. In particular, in both types of models exhibit stable cycles as a result of greater between-class competition relative to within-class competition.

[8] J.M. Cushing, Three stage semelparous Leslie models, J. Math. Biol. 59 (1) (2009) 75–104, http://dx.doi.org/10.1007/s00285-008-0208-9. [9] J.M. Cushing, A dynamic dichotomy for a system of hierarchical difference equations, J. Differ. Equ. Appl. 18 (1) (2012) 1–26. [10] J.M. Cushing, S.M. Henson, Stable bifurcations in semelparous Leslie models. J. Biol. Dyn. 6 (Suppl 2) (2012) 80–102, http://dx.doi.org/10.1080/17513758.2012. 716085. [11] J.M. Cushing, J. Li, On Ebenman’s model for the dynamics of a population with competing juveniles and adults, Bull. Math. Biol. 51 (6) (1989) 687–713. [12] N.V. Davydova, O. Diekmann, S.A. Van Gils, Year class coexistence or competitive exclusion for strict biennials? J. Math. Biol. 46 (2) (2003) 95–131, http://dx.doi. org/10.1007/s00285-002-0167-5. [13] N.V. Davydova, O. Diekmann, S.A. Van Gils, On circulant populations. I. The algebra of semelparity, Linear Algebra Appl. 398 (1–3) (2005) 185–243, http://dx. doi.org/10.1016/j.laa.2004.12.020. [14] O. Diekmann, N.V. Davydova, S.A. van Gils, On a boom and bust year class cycle, J. Differ. Equ. Appl. 11 (4–5) (2005) 327–335, http://dx.doi.org/10.1080/ 10236190412331335409. [15] O. Diekmann, S.A. van Gils, On the cyclic replicator equation and the dynamics of semelparous populations, SIAM J. Appl. Dyn. Syst. 8 (3) (2009) 1160–1189. [16] R. Dilao, A. Lakmeche, On the weak solutions of the McKendrick equation: existence of demography cycles, Math. Model. Nat. Phenom. 1 (01) (2006) 1–30. [17] S. Elaydi, An Introduction to Difference Equations, Springer Science & Business Media, 2005. [18] H.v. Foerster, Some remarks on changing populations, The Kinetics of Cellular Proliferation, (1959), pp. 382–407. [19] M.E. Gurtin, R.C. MacCamy, Non-linear age-dependent population dynamics, Arch. Ration. Mech. Anal. 54 (3) (1974) 281–300. [20] M.E. Gurtin, R.C. MacCamy, Some simple models for nonlinear age-dependent population dynamics, Math. Biosci. 43 (3–4) (1979) 199–211, http://dx.doi.org/10. 1016/0025-5564(79)90049-X. [21] M. Gyllenberg, Mathematical aspects of physiologically structured populations: the contributions of J. A. J. Metz, J. Biol. Dyn. 1 (1) (2007) 3–44, http://dx.doi.org/10. 1080/17513750601032737. [22] F. Hoppensteadt, Mathematical Theories of Populations: Deomgraphics, Genetics, and Epidemics, 20 SIAM, 1975. [23] R. Kon, Nonexistence of synchronous orbits and class coexistence in matrix population models, SIAM J. Appl. Math. 66 (2) (2005) 616–626. [24] R. Kon, Permanence induced by life-cycle resonances: the periodical cicada problem. J. Biol. Dyn. 6 (2) (2012) 855–890, http://dx.doi.org/10.1080/17513758. 2011.594098. [25] R. Kon, Y. Iwasa, Single-class orbits in nonlinear Leslie matrix models for semelparous populations, J. Math. Biol. 55 (5–6) (2007) 781–802, http://dx.doi.org/ 10.1007/s00285-007-0111-9. [26] N. MacDonald, Time lags in biological models, Lecture Notes in Biomath. 27 Springer Berlin Heidelberg, 1978. [27] A.G. McKendrick, Applications of mathematics to medical problems, Proc. Edinb.

Conflicts of interest none Funding This work was supported by the NSF grant DMS-1407564. Acknowledgments This work was completed at the University of Arizona. I would like to thank Jim Cushing for his many insightful comments on this project. References [1] A.J. Lotka, Elements of Physical Biology, Williams and Wilkins Company (1925) 435, http://dx.doi.org/10.2105/AJPH.15.9.812-b. [2] L.J.S. Allen, A density-dependent Leslie matrix model, Math. Biosci. 95 (2) (1989) 179–187, http://dx.doi.org/10.1016/0025-5564(89)90031-X. [3] H. Behncke, Periodical cicadas, J. Math. Biol. 40 (5) (2000) 413–431. [4] M.G. Bulmer, Periodical insects, Am. Nat. 111 (982) (1977) 1099–1117, http://dx. doi.org/10.2307/2678832. [5] J.M. Cushing, Bifurcation of time periodic solutions of the McKendrick equations with applications to population dynamics, Comput. Math. Appl. 9 (3) (1983) 459–478. [6] J.M. Cushing, An Introduction to Structured Population Dynamics, 71 SIAM, 1998. [7] J.M. Cushing, Nonlinear semelparous leslie models, Math. Biosci. Eng. 3 (1) (2006) 17.

10

Mathematical Biosciences 297 (2018) 1–11

A. Veprauskas

[30] R. Seydel, Practical Bifurcation and Stability Analysis, 5 Springer Science & Business Media, 2009. [31] W.O. Tschumy, Competition between juveniles and adults in age-structured populations, Theor. Popul. Biol. 21 (2) (1982) 255–268, http://dx.doi.org/10.1016/ 0040-5809(82)90017-X.

Math. Soc. 44 (1926) 98–130, http://dx.doi.org/10.1017/S0013091500034428. [28] J.A. Metz, O. Diekmann, The Dynamics of Physiologically Structured Populations, 68 Springer, 2014. [29] R. Rudnicki, R. Wieczorek, On a nonlinear age-structured model of semelparous species, Discret. Contin. Dyn. Syst. - Ser. B 19 (8) (2014) 2641–2656.

11