Bifurcations and chaos in a predator–prey model with delay and a laser-diode system with self-sustained pulsations

Abstract

Hopf bifurcations in two models, a predator–prey model with delay terms modeled by “weak generic kernel aexp(−at)” and a laser diode system, are considered. The periodic orbit immediately following the Hopf bifurcation is constructed for each system using the method of multiple scales, and its stability is analyzed. Numerical solutions reveal the existence of stable periodic attractors, attractors at infinity, as well as bounded chaotic dynamics in various cases. The dynamics exhibited by the two systems is contrasted and explained on the basis of the bifurcations occurring in each.

Introduction

Population models for single-species, two-species, and multi-species communities are of relevance in various fields of mathematical biology and mathematical ecology. Extensive reviews of various continuous models that have been investigated (such as the regular logistic equation for one species, and Lotka–Volterra systems for two species and including various additional effects as in the Kolmogorov, May, Holling, Hsu, Leslie, Caperon and other models) may be found, for example, in [1], [2]. Discrete models are considered there.

In order to incorporate various realistic physical effects that may cause at least one of the physical variables to depend on the past history of the system, it is often necessary to introduce time-delays into the governing equations. Factors that introduce time lags may include age structure of the population (influencing the birth and death rates), maturation periods (thresholds), feeding times and hunger coefficients in predator–prey interactions, reaction times, food storage times, and resource regeneration times. Models incorporating time delays in diverse biological models are extensively reviewed by MacDonald [3] and, in the context of predator–prey models, by Cushing [1].

Consider, for instance, the (modified) Lotka–Volterra two-species model

$$\text{N}\text{̇}\text{(t)=εN(1−N/κ)−αNP,}$$

$$\text{P}\text{̇}\text{(t)=−γP+βNP,}$$

where N(t) and P(t) are the prey and predator populations (numbers) respectively, ε is the birth rate of the prey, κ>0 is the carrying capacity, α is the rate of predation per predator, γ is the death rate of the predator, β is the rate of the prey’s contribution to the predator growth, and the overdot denotes a time derivative. Farkas [4] has considered the modification of (1) when a distributed (Volterra or convolution-type [5]) delay was introduced in the second equation, yielding

$$\text{N}\text{̇}\text{(t)=εN(1−N/κ)−αNP,}$$

$$\text{P}\text{̇}\text{(t)=−γP+βP∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ,}$$

where $\text{G}\text{(U)}$ is the memory function or delay kernel [3], [5]. Similarly, El-Owaidy and Ammar [6] have considered (1) modified to include delays in both equations, i.e.,

$$\text{N}\text{̇}\text{(t)=εN−αNP−}\text{εN}\text{κ}\text{∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ,}$$

$$\text{P}\text{̇}\text{(t)=−γP+βP∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ.}$$

One choice for the memory function is the Dirac delta function, which is sometimes referred to as the unit impulse function. The defining properties of this function are

$$\text{δ(t−a)=0,}\text{t≠a,}$$

$$\text{∫}{}_{\mathrm{-\infty }}{}^{\mathrm{\infty }}\text{δ(t−a)}\text{d}\text{t=1.}$$

For $\text{G}\text{(U)=δ(U)}$, one gets the case of the so-called discrete delay [1], [2], [3], [5]. While both (2) and (3) may be treated in the setting of the theory of functional differential equations [1], [3], [7], it has proven fruitful [3], [8] to employ the “linear chain” form ([3] discusses the origin of this in the theory of elasticity and in the work of Vogel [7])

$$\text{G}\text{(U)=G}{}_{\mathrm{a}}{}^{\mathrm{p}}\text{(U)≡}\text{a}{}^{\mathrm{p+1}}\text{U}{}^{\mathrm{p}}\text{p!}\text{e}{}^{\mathrm{-aU}}$$

for the memory function. Notice that in this model 1/a measures the influence of the past and as a increases this influence decreases. The functions G₁⁰(U) and G_a¹(U) (which are called the ‘weak’ and ‘strong’ kernels [1]) have been used [8] in the context of predator–prey models, such as May’s model [1].

In [4], [6], the authors used the memory function G_a⁽⁰⁾(U) to investigate the systems (2) and (3) respectively. Use of G_a⁽⁰⁾(U) reduces the integrodifferential systems (2) and (3) to the differential systems:

$$\text{N}\text{̇}\text{(t)=εN(1−N/κ)−αNP,}\text{P}\text{̇}\text{(t)=−γP+βPQ,}$$

$$\text{Q}\text{̇}\text{(t)=a(N−Q),}$$

and

$$\text{N}\text{̇}\text{(t)=εN(1−Q/κ)−αNP,}\text{P}\text{̇}\text{(t)=−γP+βPQ,}$$

$$\text{Q}\text{̇}\text{(t)=a(N−Q),}$$

respectively, where

$$\text{Q}\text{̇}\text{(t)≡∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ}$$

$$\text{=∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)a}\text{e}{}^{\mathrm{-a(t-\tau )}}\text{d}\text{τ.}$$

In particular, Farkas [4] showed that at

$$\text{a=a}{}_0\text{≡βκ−γ−ε/βκ}$$

a supercritical Hopf bifurcation [9] takes place for the system (6) and the bifurcating closed paths are asymptotically stable for a<a₀. A companion theorem was also derived, which provides simpler sufficient, but not necessary, conditions for the stability of the bifurcating closed periodic orbits. El-Owaidy and Ammar [6] showed analogously that for (7), a supercritical Hopf bifurcation takes place at a=a₀≡βκ−γ>0 (if βκ−γ>0), and the bifurcating closed periodic orbits are orbitally asymptotically stable for a<a₀.

In this paper, we will extend the analyses of Refs. [4], [5], [6] to the general predator–prey model with distributed delay in both prey and predator equations:

$$\text{N}\text{̇}\text{(t)=NF(N)−αNP−}\text{ε}\text{̃}\text{N}\text{κ}\text{∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ,}$$

$$\text{P}\text{̇}\text{(t)=−PG(P)+βP∫}{}_{\mathrm{-\infty }}{}^{\mathrm{t}}\text{N(τ)}\text{G}\text{(t−τ)}\text{d}\text{τ.}$$

For F(N)=ε=constant, G(P)=γ=constant, and $\text{ε}\text{̃}\text{=ε}$, this model would reduce to (3). For $\text{ε}\text{̃}\text{=0}$, F(N)=ε(1−N/κ) and G(P)=γ=constant, the model (9) reduces to (2). Notice that qualitative features of such general models have been considered earlier, for instance for the Kolmogorov model without delay [2] and the May model with delay [1] using G_a⁰(U) and G_a¹(U). Here, however, we will follow Refs. [4], [6], [9], [10], using G_a⁰(U) to consider the stability of the fixed points (equilibria) and the Hopf bifurcations of (9) for the general functions F(N) and G(P). This is done in Section 2. An analogous treatment of these issues for the special case of (9) with $\text{ε}\text{̃}\text{=0}$ has been considered by Roos [11]. In Section 3 we consider (9) for specific choices of F(N) and G(P) to determine the regions of phase-space where the system is dissipative (volume contracting) or dilatory (volume expanding). Section 4 considers the stability of physically relevant equilibria and Hopf bifurcation points for specific parameter values and choices of F(N) and G(P). Possible chaotic regimes are also delineated there. The systems are numerically integrated and chaotic regimes are characterized by computing power spectra, correlation functions and fractal dimensions [12].

In the remainder of this paper we shall consider Hopf bifurcations in two systems. First, we will consider (9) with $\text{G}\text{(U)=G}{}_{\mathrm{a}}{}^0\text{(U)}$, which yields the differential system

$$\text{N}\text{̇}\text{(t)=NF(N)−αNP−}\text{ε}\text{̃}\text{QN}\text{κ}\text{,}\text{P}\text{̇}\text{(t)=−PG(P)+βPQ,}$$

$$\text{Q}\text{̇}\text{(t)=a(N−Q),}$$

with Q defined by (8). This system has been investigated earlier [12] for functions F and G different from those considered here, although the primary emphasis there was on characterizing chaotic regimes using numerical diagnostics. Here, we shall carefully consider the periodic solution resulting from a Hopf bifurcation. We shall derive analytic expressions for the periodic orbit using the method of multiple scales. The stability of this periodic orbit is then considered, revealing secondary bifurcation of the orbit leading to chaos. Also, note that F(N) is chosen to incorporate the prey birth rate. Later we will consider, as a typical example, the function

$$\text{F(N)=1−δN.}$$

Similarly, the function G(P) incorporates the predator death rate and is chosen so that this rate increases with predator density P. A typical example, and the one we will use, is

$$\text{G(P)=γ(1+κP}{}^2\text{).}$$

In order to put these results in perspective, we shall also consider a semiconductor laser diode system [13], [14] with self-sustained pulsation (SSN). In normalized coordinates [14], the equations governing this system are

$$\text{d}\text{x}\text{d}\text{τ}\text{=−(x−x}{}_0\text{)−xz,}\text{d}\text{y}\text{d}\text{τ}\text{=−}\text{(y−y}{}_0\text{)}\text{τ}{}_{\mathrm{b}}\text{−γyz,}$$

$$\text{d}\text{z}\text{d}\text{τ}\text{=}\text{Γ}\text{(x+y−1)z,}$$

where τ is a normalized time, x and y are normalized carrier densities in different sections of the diode, z is a normalized photon number density, $\text{τ}{}_{\mathrm{b}}$ is a normalized carrier lifetime, γ is the ratio of the differential gains in the two sections, $\text{Γ}$ is a photon loss coefficient, and x₀ and y₀ are normalized photon pumping rates in the two sections. Physically, the parameters that may be varied most readily are the pumping rates x₀ and y₀. The other parameters, $\text{τ}{}_{\mathrm{b}}$, γ and $\text{Γ}$, are material properties of the semi-conductor and, thus, difficult to vary.

Hopf bifurcations also occur in this system leading to periodic orbits that are, once again, constructed analytically by the method of multiple scales. Stability analysis of these periodic orbits reveals that, much like the predator–prey system, they are destroyed by secondary bifurcations as the parameter is varied further past the bifurcation value. In contrast to the predator–prey system, these secondary bifurcations do not lead to chaos. This conclusion is verified numerically, showing that while the periodic orbit serves as a vehicle for SSN, which keep the system operating in a microwave oscillator mode, it does so over a narrow range of parameters.

The remainder of this paper is organized as follows. In Section 2.1 we consider the stability of the fixed point of , , and the onset of instability via a Hopf bifurcation, which may be either supercritical or subcritical. Analogous results for the laser system (13) are contained in Section 2.2. In Section 3, we derive analytical expressions for the periodic orbits resulting from this Hopf bifurcation by employing the method of multiple scales. Section 4 considers numerical solutions and discusses the results.