bifurcations in a host-parasite model with ... - Semantic Scholar
International Journal of Bifurcation and Chaos, Vol. 16, No. 11 (2006) 3291–3307 c World Scientific Publishing Company
BIFURCATIONS IN A HOST-PARASITE MODEL WITH NONLINEAR INCIDENCE WENDI WANG∗ Department of Mathematics, Southwest Normal University, Chongqing 400715, P. R. China YI LI and H. W. HETHCOTE Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA Received May 20, 2005; Revised October 26, 2005 A host-parasite model is proposed that incorporates a nonlinear incidence rate. Under the influence of multiple infectious attacks, the model admits bistable regions such that the infection dies out if initial states lie in one region, and the population and parasites coexist if initial states lie in the other region. It is also found that parasites can drive the population to extinction for suitable parameters. It is verified that the model has a saddle-node bifurcation, Hopf bifurcations and a cusp of codimension 2 or higher codimension. Stable limit cycles and unstable limit cycles are examined as the infection-reduced reproduction rate varies. It is shown that the model goes through the change of stages of infection extinction, infection persistence, infection extinction, and the extinction of both parasites and the population as the contact coefficient increases. Keywords: Nonlinear incidence; host extinction; periodic solutions; homoclinic orbits; bistable.
1. Introduction Classical epidemic models are extended in many ways to understand mechanisms of disease transmissions (see, e.g., [Hethcote, 2000; Wang & Zhao, 2004; Wang & Ruan, 2004]). An interesting topic in the study of infectious diseases is to understand how parasites regulate host populations (see [May & Anderson, 1979; De Jong et al., 1995; Grenfell et al., 1995; McCallum, 1995]). Recently, with the aim of understanding how six microparasites regulate Daphnia populations and drive the populations to extinction, Ebert et al. [2000] formulated the following epidemiological microparasite ∗
model with horizontal transmission: S+I dS = r(S + σI) 1 − − dS − λSI , dt K
(1)
dI = λIS − ( + d)I, dt where S(t) and I(t) represent the densities of uninfected (susceptible) and infected (infective) hosts at time t, respectively, r is the per capita growth rate of uninfected hosts, σ is the relative fecundity of an infected host with 0 ≤ σ ≤ 1, K is the carrying capacity of the environment for the host population, d is the parasite mortality rate, λ is the infection
Research is supported by the National Natural Science Fund of China (No. 10571143). 3291
3292
W. Wang et al.
rate coefficient, and is the infection-induced death rate. Model (1) uses the classical mass action incidence λIS and has the usual asymptotic behavior. If there is no endemic equilibrium, the infection free equilibrium is globally stable. If it does have an infected equilibrium, this steady state is globally stable. Thus this model fails to explain the observed rich outcomes that depend on parameter values and initial population levels. To explain why the host population can be driven to extinction, a carefully designed stochastic simulation of model (1) was conducted in [Ebert et al., 2000]. The simulations indicated that extinction of the host is possible in some parameter regions. This means that stochastic factors are possible causes of the extinction of the host population. The mass-action principle stems from the theory of the kinetics of chemical reactions. It implies that individuals, like molecules, would be mixed homogeneously in space and every encounter would have the same occurrence probability. This may be applicable at low host population densities. For large densities, the number of new encounters made by a single individual per unit time cannot increase linearly with population sizes, either because overcrowding reduces the movements of individuals, or because a related increase of occupied habitat prevents its total occupation by individuals. In this case, a standard incidence λSI/(S + I) would be a suitable option (see [Hethcote, 2000] and [Deredec et al., 2003]). Diekmann and Heesterbeek [2000] gave a more general incidence C(N )SI/N with N = S +I, where C(N ) is the number of individuals that are engaged in contacts at any time. This incidence is suitable for any population size, and reduced to the mass action law if C(N ) is linearly dependent on N at low population densities, or reduced to the standard incidence if C(N ) takes a constant value at high population densities. However, the standard incidence may also apply to low population densities due to adaptive behavior of populations. It is observed that larger populations may occupy larger areas [Begon et al., 2002]. This phenomenon seems more common for wild populations. Based upon these observations, [De Jong et al., 1995] and [Bouma et al., 1995] suggested that each individual occupies a characteristic area such that the area occupied by a population is directly proportional to the size of that population. In this case, the incidence at low host population size may also be represented or approximated by the standard incidence.
Hwang and Kuang [2003] changed the mass action incidence in (1) to the standard incidence λSI/(S + I). Significantly, this revised parasite-host model can exhibit the observed parasite-induced host extinction. This means that the extinction of the host population may be due to deterministic factors, instead of stochastic factors. Since the mass action incidence implies that individuals of the population occupy a fixed area and the standard incidence implies that the occupied area is proportional to its population size, alternately, the extinction of the host population may be from the behavior changes of the population from the movement in a fixed area to the movement in an area that is proportional to its population size. However, periodic fluctuations of the population and outcomes that are dependent on initial population levels are not predicted in the model with the standard incidence, as pointed out in their paper [Hwang & Kuang, 2003]. Note that the standard incidence can be further written as β0 ωSI/(S + I) where ω is the encounter rate and β0 is the probability that a susceptible individual becomes infected in an encounter by an infectious individual. The standard incidence implies that the infection probability is invariant for all sizes of infectious individuals and susceptibility of the individuals. However, there is evidence that infections are limited by host physiological resistance and behavioral defense of the immune system [Nouhuys et al., 2004]. In these cases, a susceptible is more vulnerable under multiple attacks such that multiple attacks by infectious members may increase the infection probability of the susceptible. Thus, it is reasonable to assume that β0 is an increasing function of the number of infectious attacks ωI/(S + I) per unit time. To explore the influence of this infection probability as a function of multiple attacks on the dynamics of hostparasite evolutions, in this paper we consider only the simplest case that the infection probability is proportional to the number of infectious attacks ωI/(S + I). Then we have β0 = β1 ωI/(S + I), where β1 is a proportional constant. As a consequence, the incidence becomes βSI 2/(S + I)2 with β = β1 ω 2 . Indeed, this incidence is a special case of the infection force βG(I/N ) of proportional mixing which was proposed and used by Nold [1980], Derrick and van den Driessche [1993], Castillo-Chavez and Yakubu [2001] and Derrick and van den Driessche [2003]. If a population size is constant, this type of nonlinear incidences was studied by [Liu et al.,
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
1986, Liu et al., 1987; Lizana & Rivero, 1996; Ruan & Wang, 2003]. Based upon the incidence and the modeling approach for the endemic models with the logistic demographic structure in [Gao & Hethcote, 1992], we formulate the following SI endemic model: χrN dN = b− [S + σI] dt K (1 − χ)rN N − I − d+ K N N − ( + b(1 − σ))I =r 1− K χrIN K χrN = b− [S + σI] K I2 (1 − χ)rN S − βS 2 , − d+ K N I2 (1 − χ)rN = βS 2 − + d + I, N K + (1 − σ)
dS dt
dI dt
(2)
where N = S + I is the total host population density, b − χrN/K is the per capita birth rate for susceptibles, and σ is the infection-induced reduction in the birth rate, so that [b − χrN/K]σ is the per capita birth rate for infectives. Here d+(1−χ)rN/K is the natural per capita death rate, r = b − d > 0 is the per capita net growth rate, K is the environmental carrying capacity, χ is the convex combination constant with 0 ≤ χ ≤ 1, and is the per capita infection-related death rate. Note that the population dynamics are governed by the logistic differential equation, when there is no infection (I = 0). For 0 < χ < 1, the per capita birth rate b − χrN/K is a decreasing function of N and the natural per capita death rate d + (1 − χ)rN/K is an increasing function of N , which is consistent with densitydependent population dynamics. When χ = 1, the density dependence affects the per capita birth rate, but does not affect the per capita death rate. In order to get an SI epidemic model comparable to those considered by Elbert et al. [2000], Hwang and Kuang [2003], we set χ = 1, so that the model (2) becomes rN I2 dS = b− (S + σI) − dS − βS 2 , dt K N (3) 2 I dI = βS 2 − ( + d)I. dt N
3293
This formulation of an SI endemic model with the density-dependent growth is slightly different from that in model (1), in which dN N = r[N − (1 − σ)I] 1 − − dN − I. dt K Thus, in the absence of the disease, model (1) implies rN dN = r−d− N dt K so that N → (r − d)K/r as t → ∞. The advantage of the model (3) is that the original meanings of parameters r and K in the logistic differential equation are maintained, in the sense that r is the intrinsic growth rate of the population and N → K as t → ∞ in the absence of the disease. For biological reasons, we will restrict out attention on the region {(S, I) : 0 ≤ S + I ≤ K}. Note that (3) has the singularity at (0, 0). We make a time scale change dt = N 2 dτ such that (3) is equivalent to the following system in the interior of the first quadrant R2 : dS = (r + 2bσ − β)SI 2 + (2r + bσ)S 2 I + rS 3 dt + bσI 3 −
rS 4 rσI 4 r(3 + σ)S 3 I − − K K K
r(1 + σ)S 2I 2 r(3σ + 1)SI 3 − , −3 K K
(4)
dI = I[βSI − ( + d)(S + I)2 ], dt where t represents τ for the convenience of notations. To be concise in analysis, let us scale (4) by x = rS/(βK), y = rI/(βK) and ξ = β 3 K 2 t/r 2 . If t is used to represent ξ, we obtain dx = ax3 + (2a + σh)x2 y + (a + 2σh − 1)xy 2 dt + σhy 3 − x4 − (3 + σ)x3y − 3(1 + σ)x2 y 2 − (3σ + 1)xy 3 − σy 4 ,
(5)
dy = y(xy − c(x + y)2 ), dt where a=
r , β
h=
d+r , β
c=
d+ . β
In this paper, we assume that r, β, d and are positive constants, and σ is a constant in [0, 1]. Hence, we
3294
W. Wang et al.
have h > a and a + c − hσ > 0. The feasible region of (5) is {(x, y) : x ≥ 0, y ≥ 0, 0 ≤ x + y ≤ a}. The organization of this paper is as follows. In the next section, we analyze dynamical behavior of (5) around (0, 0) and give conditions under which the infection cannot drive the population to extinction, and conditions under which the infection drives the population to extinction. In Sec. 3, we present the analysis of bifurcation structures including Hopf bifurcations and higher codimension bifurcations, and show the existence of limit cycles. We conclude with a discussion of the results in Sec. 4.
2. Extinctions In this section, we analyze the existence and stability of boundary equilibria of (5) and present conditions under which the infection drives the population to extinction and conditions under which the infection cannot drive the population to extinction. First, (a, 0) is an infection free equilibrium. It is always asymptotically stable. This means that the outcome of infection is unsuccessful if the invasion of parasites is small. To be endemic, the intensity of invasion must be higher than a threshold. Secondly, E0 = (0, 0) is a trivial equilibrium. This corresponds to the extinction of both infected hosts and the host population. Since E0 is nonhyperbolic, we introduce the polar coordinates x = R cos θ, y = R sin θ to obtain dR = R3 (V (θ) + o(1)), dt
(6)
dθ = −R2 sin θ(cos θ + sin θ)(U (θ) + o(1)), dt where U (θ) = (a + c)cos2 θ + (a + c + hσ − 1)cos θ sin θ + hσ sin2 θ, V (θ) = a(cos(θ))2 (sin(θ) + cos(θ))2 + hσ sin(θ)cos(θ)(sin(θ) + cos(θ))2 − c(sin(θ))2 (sin(θ) + cos(θ))2 + (sin(θ))2 cos(θ)(− cos(θ) + sin(θ)). Note that the first quadrant R2 is positively invariant. By the qualitative theory of differential equations [Zhang et al., 1991] (see also [Xiao & Ruan, 2001; Berezovskaya et al., 2001] for population models), an orbit of (5), if it tends to the origin, must tend to it along a characteristic direction. θ = 0 is a trivial characteristic direction, where infected hosts die out and the uninfected
host population tends to its carrying capacity K. We now look for other characteristic directions in (0, π/2). If z = tan θ, by the definition of U we see that z satisfies F1 (z) := hσz 2 + (a + c + hσ − 1)z + a + c = 0. (7) Thus, there is no characteristic direction in (0, π/2) if a + c + hσ − 1 > 0 or
a + c + hσ − 1 < 0, (a + c + hσ − 1)2 < 4hσ(a + c).
(8)
(9)
If (8) or (9) is satisfied, since V (0) = a > 0 and U (0) = a + c, by Theorems 3.5 and 3.7 of [Zhang et al., 1991], we have the following theorem: Theorem 1. If (8) or (9) is satisfied, then there is
an 1 > 0 and an R1 > 0 such that (5) admits a unique orbit in {(θ, R) : 0 ≤ θ < 1 , 0 < R < R1 } that tends to (0, 0) along θ = 0 as t → −∞. Biological implication of Theorem 1 is that the population and parasites cannot die out together because (0, 0) is a repeller. Note that (7) has two positive solutions z1 and z2 with z1 < z2 if a + c + hσ − 1 < 0, (10) (a + c + hσ − 1)2 > 4hσ(a + c). This means that we have two other characteristic directions 0 < θ1 < θ2 except for θ = 0. In order to find dynamical behaviors of (5) along these two directions, we apply the Briot–Bouquet transformation: x = x,
y = ux,
dτ = x2 dt
to transform (5) into dx = x[a − x + (2a + hσ)u − (3 + σ)xu dt + (a + 2hσ − 1)u2 − 3(1 + σ)xu2 + hσu3 − (3σ + 1)xu3 − σxu4 ], du = −u(u + 1)[a + c − x dt
(11)
+ (a + c + hσ − 1)u − (2 + σ)xu + hσu2 − (1 + 2σ)xu2 − σxu3 ], where t is used to represent τ . We consider (11) in the region x ≥ 0, u ≥ 0 due to the biological
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
background. Under the condition (10), we see that (11) admits equilibria (0, 0), (0, z1 ) and (0, z2 ). It is easy to see that (0, 0) is a saddle. To analyze equilibrium (0, z1 ), we make the change of variables x1 = x, x2 = u − z1 to transform (11) to dx1 = [a + (2a + hσ)z1 + (−1 + a + 2hσ)z12 dt + hσz13 ]x1 + o((x1 , x2 )), dx2 = z1 (1 + z1 )3 (σz1 + 1)x1 + [−a − c dt + (−4a − 2hσ − 4c + 2)z1 (−6hσ + 3 − 3a − 3c)z12 − 4hσz13 ]x2 + o((x1 , x2 )),
(12)
Using (7), (12) can be simplified to dx1 = g1 (z1 )x1 + o((x1 , x2 )), dt dx2 = z1 (1 + z1 )3 (σz1 + 1)x1 dt + g2 (z1 )x2 + o((x1 , x2 )),
Further, (0, z1 ) of (11) is a saddle if hσ < h3 σ < c(1 + (a/(c − 1))) and F2 (h) > 0. Note that A0 := a + c − hσ > 0. Condition (i) means
Proof.
A1 := −cσh − c + ca + c2 + hσ ≥ 0. Hence, g1 (z1 ) = A1 z1 + A0 > 0. If (ii) holds, then A1 < 0. Since h ≥ h3 , we have F0 (h) < 0. By direct calculations, we see that F0 (h) has the same sign as A0 a + c + hσ − 1 . − A1 2hσ Hence, A0 > 0. g1 (z1 ) = A1 z1 − −A1 When (iii) holds, we have a + c + hσ − 1 A0 . > A1 2hσ
(13)
where g1 (z1 ) = c(a − hσ + c) + (−cσh − c + ca + c2 + hσ)z1 , g2 (z1 ) = −(a + c)(a − hσ + c − 1) + (hσ(c + a + 1) − (a + c − 1)2 )z1 . The stability of (0, z1 ) is determined by the signs of gi (z1 ). Let us consider (a + c + hσ − 1)2 − 4 hσ(a + c) = 0.
3295
(14)
This equation in h has two solutions 0 < h1 < h2 . It is easy to verify that (10) implies 0 < h < h1 . Next, we consider F0 (h) := h2 σ 2 (c + 1) − hσ(−a + 2c2 − c + 2ac + 1) + c(a + c − 1)2 = 0. (15) This equation in h has two solutions 0 < h3 < h4 with h1 < h4 < h2 . Set F2 (h) = ch2 σ 2 − hσ(2c2 − a + 2ac) + c(c2 + 2ac + a2 − a).
Note that A0 hσ = F1 −A1 (−σhc − c + ac + c2 + hσ)2 × F2 (h) < 0. Thus, A0 /(−A1 ) > z1 . Consequently, we have g1 (z1 ) > 0. Now, we consider the sign of g2 (z1 ). Since a + c + hσ < 1, we have B0 := −(a + c)(a − hσ + c − 1) > 0. Set B1 := hσ(c + a + 1) − (a + c − 1)2 . Using the fact that h < h1 , we obtain hσ < (a + c − 1)2 /(a + c + 1), which means B1 < 0. Then by similar arguments to above, we obtain a + c + hσ − 1 B0 , F2 (h) < 0, > B1 2hσ where h < h1 is used. Thus, B0 /(−B1 ) > z1 . It follows that B0 > 0. g2 (z1 ) = B1 z1 − −B1
Theorem 2. Let (10) hold. Then (0, z1 ) of (11) is an unstable node if one of the following conditions is satisfied:
In summary, we have verified that g1 (z1 ) > 0 and g2 (z1 ) > 0 if any of (i)–(iii) is satisfied. Thus, (0, z1 ) is an unstable node in this case. The conclusion that (0, z1 ) is a saddle can be proved in a similar way.
(i) c(1 + (a/(c − 1))) < hσ; (ii) h3 σ ≤ hσ < c(1 + (a/(c − 1))); (iii) hσ < h3 σ < c(1 + (a/(c − 1))) and F2 (h) < 0.
By similar discussions to those in the proof of Theorem 2, we can obtain the following results for the stability of (0, z2 ).
3296
W. Wang et al.
0.4
0.35
0.3
I
0.25
0.2
0.15
0.1
0.05
0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
S
0.1
0.09
0.08
0.07
I
0.06
0.05
0.04
0.03
0.02
0.01
0 0
0.05
0.1
0.15
S
Fig. 1.
The figures of two typical ways of extinction of both host population and parasites.
Theorem 3. Let (10) hold. Then (0, z2 ) of (11) is a saddle if one of the following conditions is satisfied:
(i) c(1 + (a/(c − 1))) < hσ; (ii) h3 σ < hσ < c(1 + (a/(c − 1))) and F2 (h) > 0. Further, (0, z2 ) of (11) is a stable node if one of the following holds: (1) h3 σ < hσ < c(1 + (a/(c − 1))) and F2 (h) < 0; (2) hσ < h3 σ < c(1 + (a/(c − 1))). From Theorems 2 and 3, there are three types of dynamical behaviors around (0, 0) for (5) when the characteristic directions θ1 and θ2 exist. First, if (0, z1 ) is an unstable node and (0, z2 ) is a saddle, then (0, 0) of (5) is a repeller so that the host population and parasites cannot die out together.
Secondly, if (0, z1 ) is a saddle and (0, z2 ) is a stable node, the region in the first quadrant around (0, 0) is split into two parts, one of which is an attracting parabolic sector and the other is a hyperbolic sector. The host population and parasites go to extinction in the attracting parabolic sector. Lastly, if (0, z1 ) is an unstable node and (0, z2 ) is a stable node, the region in the first quadrant around (0, 0) is split into two parts, one of which is a repelling parabolic sector and the other is an elliptic sector. The host population and parasites go to extinction in the elliptic sector (see Fig. 1).
3. Bifurcations 3.1. Positive equilibria We now consider positive equilibria of (5). Roughly speaking, the existence of these equilibria is
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
necessary for the coexistence of the population and parasites. First, we perform a topological transformation to (5) so that positive equilibria can be easily found. Set u = x + y, v = y/(x + y) and dτ = x2 dt. Substituting them into (5) and then replacing u by x, v by y, τ by t for the simplicity of notations, we obtain dx = x(a − x + (hσ − c − a)y + (1 − σ)xy), dt dy = y(−c − a + x + (1 + c + a − hσ)y dt − (1 − σ)xy − y 2 ).
(16)
Evidently, we should consider this system in the region D = {(x, y) : 0 < x < a, 0 < y < 1}. The isoclines of (16) in D are given by x= x=
a − (a − hσ + c)y = f1 (y), 1 − (1 − σ) y
a + c − (1 + c + a − hσ)y + 1 − (1 − σ)y
y2
= f2 (y).
(17) (18)
(19)
If hσ < c, we have a/(a + c − hσ) < 1. Thus, it suffices to consider solutions of (19) in (0, a/(a + c − hσ)]. Note that f3 (0) = c > 0 and F2 (h) a = . f3 a + c − hσ (a + c − hσ)2
x(−σy + y − 1) J := y(σy − y + 1)
If c ≥ hσ + a, then a/(a + c − hσ) ≤ 1/2. Note that f3 is a decreasing function on (0, a/(a + c − hσ)). Consequently, (19) has a unique solution on (0, a/(a + c − hσ)) if F2 (h) < 0, and has no solution on (0, a/(a + c − hσ)) if F2 (h) > 0. Similarly, if hσ < c < hσ + a, then (19) has a unique solution on (0, a/(a + c − hσ)) if c < 1/4 and F2 (h) < 0, has two solutions on (0, a/(a + c − hσ)) if c < 1/4 and F2 (h) > 0, and has no solution on (0, a/(a+c−hσ)) if c > 1/4. For the case where c ≤ hσ, since f3 (1) = c > 0, it follows that (19) has two positive solutions on (0, 1) if c < 1/4, has a unique positive solution on (0, 1) if c = 1/4 and has no solution on (0, 1) if c > 1/4. In summary, we have the following results for positive equilibria of (16). Theorem 4. (16) has a unique positive equilibrium
E ∗ if one of the following conditions holds:
To obtain a positive equilibrium in D, we should have y < a/(a + c − hσ). By (17) and (18), we obtain f3 (y) := y 2 − y + c = 0.
3297
(i) c ≥ hσ + a and F2 (h) < 0; (ii) hσ < c < min{hσ + a, 1/4} and F2 (h) < 0; (iii) c = 1/4 < hσ. Further, (16) has two positive equilibria E1 = (x1 , y1 ) and E2 = (x2 , y2 ) where y1 < 1/2, y2 > 1/2 and xi = f1 (yi ), i = 1, 2 if one of the following conditions is satisfied: (1) hσ ≤ c < min{hσ + a, 1/4} and F2 (h) > 0; (2) c < min{hσ, 1/4}. Now, we consider the stability of positive equilibria. The Jacobian matrix at a positive equilibrium (x, y) is
−x(xσ − x + c + a − hσ) . y(1 − 2y + xσ − x + c + a − hσ)
Its characteristic equation is λ2 + (σyh − ay − yc + x − y + 2y 2 )λ + yx(2y − 1)(σy − y + 1) = 0. Thus, det(J) < 0 if y < 1/2 and det(J) > 0 if y > 1/2. Consequently, we can state. Theorem 5. If (16) has a unique endemic equilib-
rium E ∗ and (iii) in Theorem 4 does not hold, this endemic equilibrium is a saddle. Further, if (16) has two endemic equilibria, then E1 is a saddle, E2 is a node or focus. One implication of Theorems 4 and 5 is the existence of a saddle node bifurcation. For example, if
we fix d = 0.2, = 0.1, r = 0.16, σ = 0.6, then a = 0.16/β, h = 0.36/β, c = 0.3/β. In this case, there are two positive equilibria when 1.2 < β < 1.3289, one of which is a saddle and the other one is a node or a focus (see Fig. 2). To indicate another way of occurrence of positive equilibria, we set d = 0.2, = 0.1, r = 0.16, β = 1.22. Then a = 0.1311, h = 0.2951, c = 0.2459. It is easy to see that (16) does not have a positive equilibrium if 0 ≤ σ < 0.25, has a unique positive equilibrium if 0.25 < σ < 0.4898. When σ increases from 0.4898, (0, 0) becomes a repeller and the second positive equilibrium emerges.
W. Wang et al.
3298
If a1 < 0 and a0 ≤ 0, then the trace is negative. Therefore, E2 is stable. If a1 < 0 and a0 > 0, we have a0 tr 1 = a1 y2 − . (20) −a1
0.07 0.06 0.05 0.04
x
Evidently, tr 1 < 0 if 0.03
2a0 + a1 < 0. 0.02 0.01 0 −0.01 1.18
1.2
Fig. 2.
1.22
1.24
1.26 β
1.28
1.3
1.32
1.34
The graph of a saddle node bifurcation.
We now determine the stability of E2 = (x2 , y2 ). Set H = a + c − hσ and a1 = H(1 + σ) − 2c − σ + 2cσ a0 = H(−1 + c − cσ) − hσ + 2c + cσ, r1 = (c2 (1 − σ)2 + σ(4c − 1))H 2 + H(hσ(−1 + 2c)(−1 + σ) + σ − 4cσ − 2c2 + 2c2σ) + h2σ 2 − hσ(2c + 4cσ − σ) + c(−2σ − cσ 2 + 10cσ + 4c2 (1 − σ)2 ). Theorem 6. Assume that (16) has two endemic equilibria. Then E2 is asymptotically stable if any one of the following conditions is satisfied
(1) a1 < 0 and 2a0 + a1 < 0; (2) a1 < 0, 2a0 + a1 > 0, r1 < 0; (3) a1 > 0, 2a0 + a1 < 0 and r1 > 0. Further, E2 is unstable if any one of the following conditions is satisfied (i) a1 > 0 and 2a0 + a1 > 0; (ii) a1 > 0, 2a0 + a1 < 0 and r1 < 0; (iii) a1 < 0, 2a0 + a1 > 0 and r1 > 0. Note that the determinant of the Jacobian matrix J is positive at E2 . It suffices to consider the sign of the trace of J. Note that the sign of the trace is determined by Proof.
(−x + y − σhy2 + ay2 +
cy2 − 2y22 ).
Using (19) and (17), we see that its sign is determined by tr 1 := a1 y2 + a0 .
(21)
Hence, E2 is stable if condition (1) holds. By direct calculations, we see that the sign of f3 (−a0 /a1 ) is determined by r1 . If a1 < 0, a0 > 0 and 2a0 + a1 > 0, which implies −a0 /a1 > 1/2, it follows from r1 < 0 that −a0 /a1 < y2 , which leads to tr 1 < 0. Hence, E2 is stable if condition (2) holds. Similarly, E2 is stable if condition (3) is satisfied. The proofs that E2 is unstable when one of the conditions (i)–(iii) holds are similar to the above discussions. At this stage, we are able to indicate certain global behaviors of (5). First, we consider the case where (5) does not have a positive equilibrium. Then there are two types of dynamical behaviors. If (0, 0) is a repeller, the infection-free equilibrium is globally stable, which is stated in Theorems 7 and 8. If there is an attracting sector, two split regions occur where the host population and parasites die out in one region and the infection dies out in the other region. This case is described in Theorem 9. Theorem 7. Let (8) or (9) hold. Then the infection-
free equilibrium (a, 0) is globally stable if one of the following conditions is satisfied: (i) c > hσ + a and F2 (h) > 0 holds; (ii) c < hσ + a and c > 1/4. Theorem 1 implies that (0, 0) repels positive orbits if (8) or (9) holds. By the arguments preceding Theorem 4, we see that (5) does not have a positive equilibrium if (i) or (ii) holds. It is easy to verify that positive solutions of (5) are ultimately bounded. It follows that positive solutions of (5) tend to the infection-free equilibrium. Since it is always stable, the global stability follows.
Proof.
Similarly, we can obtain the following results: Theorem 8. Let (10) hold. Assume that one of the
conditions (i)–(iii) in Theorem 2 is satisfied and that one of the conditions (i) and (ii) in Theorem 3 is satisfied. Then the infection-free equilibrium (a, 0) is globally stable if one of the following conditions
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
is satisfied:
all positive orbits of (5) initiating from the other region tend to (0, 0) along the characteristic direction θ2 as t tends to infinity.
(i) c > hσ + a and F2 (h) > 0 holds; (ii) c < hσ + a and c > 1/4. Theorem 9. Let (10) hold. Assume further that
one of the conditions (1) and (2) in Theorem 3 is satisfied. If one of the following conditions:
Example 11. Let us consider the differential equa-
tions in (3). r(S + I) dS = (S + σI) b − − dS dt K 2 I S, −β S+I 2 I dI =β S − ( + d)I. dt S+I
(i) c > hσ + a and F2 (h) > 0 holds, (ii) c < hσ + a and c > 1/4, is satisfied, then there is a positive orbit in (5) which tends to (0, 0) along the characteristic direction θ1 as time tends to ∞ or −∞, and which splits the interior of (5) into two regions such that all positive orbits of (5) initiating from one region tend to (a, 0) as t tend to infinity, and all positive orbits of (5) initiating from the other region tend to (0, 0) along the characteristic direction θ2 as t tends to infinity. We have similar results if there is a unique positive equilibrium which is nondegenerate. Theorem 10. Let (10) hold. Assume further that
one of the conditions (1) and (2) in Theorem 3 is satisfied. If one of the following conditions: (i) c ≥ hσ + a and F2 (h) < 0; (ii) hσ < c < min{hσ + a, 1/4} and F2 (h) < 0; is satisfied, then the stable manifolds of E ∗ splits the interior of feasible region of (5) into two regions such that all positive orbits of (5) initiating from one region tend to (a, 0) as t tend to infinity, and
0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 3.
0.2
0.4
0.6
0.8
(22)
We fix r = 0.5, d = 0.4, σ = 0.1, = 0.5. After we transform (22) into (5), we obtain a = 0.5/β, h = 0.9/β and c = 0.9/β. When 0 < β < 2.1999, Theorem 7 is satisfied. Thus, the infectionfree equilibrium is globally stable. When 2.1999 < β < 3.8136, Theorem 9 is satisfied where there is no endemic equilibrium. When β > 3.8136, Theorem 10 is satisfied where there is a unique endemic equilibrium. Thus, in the last two cases, there are bistable regions and the outcome of evolution of the population and parasites depends on their initial conditions (see Figs. 3 and 4. This figure and others in the following are produced by the package PPLANE6 [Polking, 2003] unless stated otherwise). The infection dies out and the population tends to its carrying capacity in the lower region, while both the population and parasites die out in the upper region.
1
0
3299
1 S
1.2
1.4
1.6
1.8
2
The graph of bistable regions where r = 0.5, d = 0.4, σ = 0.1, = 0.5, K = 2 and β = 2.5.
3300
W. Wang et al. 1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 4. The stable manifolds of E ∗ split the feasible region into two parts. The lower part is a basin of the infection-free equilibrium and the upper part is a basin of (0, 0), where r = 0.5, d = 0.4, σ = 0.1, = 0.5, K = 2 and β = 4.5.
We should address that there is an elliptic sector in the case where Theorem 10 is satisfied (see Fig. 4). When two positive equilibria occur, global structure of (5) is more complicated because periodic solutions and homoclinic orbits may occur. To obtain information of periodic solutions of (16), we consider its Hopf bifurcations.
Theorem 12. Let condition (1) or condition (2) in Theorem 4 hold and r1 = 0. Assume that a1 < 0, 2a0 + a1 > 0 or a1 > 0, 2a0 + a1 < 0. Then there is a stable periodic orbit in (16) near the Hopf bifurcation point when ξ > 0, and there is an unstable periodic orbit in (16) near the Hopf bifurcation point when ξ < 0.
In order to determine the direction of a Hopf bifurcation, let us make the transformation of X = x − x2 , Y = y − y2 . Then (16) becomes
Proof.
3.2. Hopf bifurcation A Hopf bifurcation may occur when c < min{hσ + a, 1/4} and r1 = 0. In order to know the direction of a Hopf bifurcation, we set ξ := ((a1 + 4a0 )(1 − σ) + 2a1 )a1 H − + (a21 + 4a1 a0 + 6a20 )(1 − σ).
a21
dX = a11 X + a12 Y + g1 (X, Y ), dt
(23)
dY = a21 X + a22 Y + g2 (X, Y ), dt
where a11 = a12 =
a0 (a0 (1 − σ) + a1 )(Ha1 + a1 + 2a0 ) , a31
(Ha 21 + (1 − σ)(a1 a0 + a1 a0 H + 2a20 ))(Ha 1 + a1 + 2a0 )a0 , a41
a21 = −
a0 (a0 (1 − σ) + a1 ) , a21
a22 = −a11 and gi , i = 1, 2, are terms of higher order. Now, let Z = a11 X + a12 Y to obtain Z − a11X dZ Z − a11X Z − a11X dX = Z + g1 X, = −a2 X + a11 g1 X, , + a12 g2 X, , dt a12 dt a12 a12 where a2 = −
(a1 + 2a0 )(a0 (1 − σ) + a1 )(Ha1 + a1 + 2a0 )a20 . a51
(24)
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
3301
√ Set u = −X, v = Z/ a2 . Then (24) becomes √ du = − a2 v + G1 (u, v), dt where
Set
dv √ = a2 u + G2 (u, v), dt
(25)
√ v a2 + a11 u , G1 (u, v) = −g1 −u, a12 √ √ v a2 + a11 u v a2 + a11 u a11 g1 −u, + a12 g2 −u, a12 a12 . G2 (u, v) = √ a2 2 2 ∂ 3 G2 ∂ 3 G1 ∂ 3 G2 1 ∂ G1 ∂ G1 ∂ 2 G1 1 ∂ 3 G1 + + 2 + + + √ Φ(u, v) = 16 ∂u3 ∂u∂v 2 ∂u ∂v ∂v 3 16 a2 ∂u∂v ∂u2 ∂v 2 ∂ 2 G1 ∂ 2 G2 ∂ 2 G1 ∂ 2 G2 ∂ 2 G2 ∂ 2 G2 ∂ 2 G2 + + − . − ∂u∂v ∂u2 ∂v 2 ∂u2 ∂u2 ∂v 2 ∂v 2
By the results in [Guckenheimer & Holmes, 1996] or [Perko, 1996], the direction of the Hopf bifurcation is determined by the sign of Φ(0, 0). The expression of Φ(0, 0) is very complicated. However, using the facts that a11 + a22 = 0 and y22 − y2 + c = 0, we have c=−
a0 (a1 + a0 ) , a21
a=−
a0 ((a1 a0 H + a1 a0 + 2a20 )(1 − σ) + 2Ha21 + 2a0 a1 + a21 ) . a31
As a consequence, we see that the sign of Φ(0, 0) is determined by ξ where ξ : = ((a1 + 4a0 )(1 − σ) + 2a1 )a1H − a21 + (a21 + 4a1 a0 + 6a20 )(1 − σ). The conclusion of Theorem 12 now follows from [Guckenheimer & Holmes, 1996] or [Perko, 1996]. Example 13. Let us fix d = 0.2, = 0.1, r =
0.16, β = 1.22. Then a = 0.1311, h = 0.2951, c = 0.2459, H = 0.3771 − 0.2951σ. If we vary σ from 0 to 1, we have a1 = −0.2951(σ + 1.0865)(σ + 0.3579) < 0, a0 = 0.0726σ 2 + 0.0806σ + 0.2075 > 0. Further, we have 2a0 + a1 > 0 if 0 ≤ σ < 0.7845 and 2a0 + a1 < 0 if 0.7845 < σ ≤ 1. Then it is easy to see that Ei , i = 1, 2, exist when 0.4898 < σ ≤ 1. Note that r1 = 0.0053(σ + 2.8579)(σ + 2.3494) × (σ − 0.6473)(σ − 0.9822). Thus, E2 is stable when 0.6473 < σ ≤ 1 and is unstable when 0.4898 < σ < 0.6473. This means that a Hopf bifurcation may occur when σ passes through 0.6473. Since ξ1 = −0.125 when σ = 0.6473, it follows from Theorem 12 that (16)
admits an unstable limit cycle when σ increases from 0.6473. If 0 ≤ σ < 0.4898, the model has no positive equilibrium or has one positive saddle point, and the dynamical behavior is similar to those indicated in Figs. 3 and 4. If 0.4898 < σ < 0.6473, computer simulations show that all positive orbits except two endemic equilibria and the stable manifolds of E1 tend to the infection-free equilibrium as t tends to ∞ (see Fig. 5). This means that the infection dies out although two endemic equilibria occur. When σ increases from 0.6473, we have an unstable limit cycle inside which all orbits approach E2 as t tends to ∞ and outside which all orbits tend to the infection-free equilibrium as t tends to ∞ (see Fig. 6). Further, as σ increases, the limit cycle expands and a homoclinic orbit occurs when σ = 0.667938 (see Fig. 7). Then, if we further increase σ, the homoclinic orbit is broken, and the stable manifolds of E1 split the feasible region into two parts. If an initial position lies in the lower part, the orbit tends to the infection-free equilibrium. If an initial position lies in the upper part, the orbit tends to the endemic equilibrium E2 (see Fig. 8). Hence, the infection can drive the population to extinction if the reduced reproduction coefficient σ is small, the infection dies out if the coefficient is medium, and the population coexists with parasites if the coefficient is close to 1. An unstable limit cycle
3302
W. Wang et al. 1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 5. E1 is a saddle and E2 is unstable. All the positive orbits except E1 , E2 and the stable manifolds of E1 tend to the infection-free equilibrium, where d = 0.2, = 0.1, r = 0.16, β = 1.22, σ = 0.55 and K = 2.
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
Fig. 6.
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
An unstable limit cycle occurs around E2 where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2 and σ = 0.66.
and a homoclinic orbit are the boundaries of the coexistence regions. Now, we consider two important cases to see the effects of the infection-related death and the infection-reduced birth on the dynamical behavior of (16). First, we suppose σ = 1. This means that parasites are harmless for the reproduction of the population. In this case, the boundary equilibria are the same as those discussed in last section, but positive equilibria and their Hopf bifurcations can be simplified. Indeed, we have a1 = −2H + 1, a0 = H + h − 3c, 2a0 + a1 = 2h − 6c + 1 > 0. By
Theorem 6, we see that a Hopf bifurcation can only occur when H > 1/2 and (4c − 1)H 2 + (1 − 4c)H − 2c − 6ch + h2 + 9c2 + h = 0.
(26)
Furthermore, we have ξ = a21 (2H − 1) > 0. Therefore, we can state Theorem 14. Let E2 exist. If σ = 1 and (26) holds,
then (16) has an unstable limit cycle near the bifurcation point.
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
3303
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
Fig. 7.
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
A homoclinic orbit occurs when σ = 0.667938 where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2.
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 8. The stable manifolds of E1 split the feasible region into two regions one of which is the basin of E2 and the other one is the basin of the infection-free equilibrium, where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2 and σ = 0.72.
Now, we consider the case where σ = 0. This case means that infected hosts do not have fecundity. Then a1 = −2c + H, a0 = 2c − H + cH, 2a0 + a1 = 2c − H + 2cH. By Theorem 6, we see that a Hopf bifurcation can only occur when H > 2c/(1 − 2c) or H < 2c where c < 1/4 is used, and 2
H − 2H + 4c = 0.
(27)
By direct calculations, we obtain ξ=− Hence, we have
H 5 (H − 2) . 8
(28)
Theorem 15. Assume that E2 exists and σ = 0.
Then (16) has an unstable limit cycle near the Hopf bifurcation point if H > 2 and has a stable limit cycle near the Hopf bifurcation point if 2c/(1 − 2c) < H < 2 or H < 2c. Example 16. Let us fix r = 1, d = 0.1, σ =
0, = 0.1, β = 0.9 in (3). Then we have a = 1.1111, h = 1.2222, c = 0.2222, H = 1.3333. It is easy to see that condition (1) in Theorem 4 holds. Since 0.5714 = 2c/(1 − 2c) < H < 2, it follows that (3) admits a stable limit cycle when β is close to 0.9 (see Fig. 9). By computer simulations, we
3304
W. Wang et al. 1.5
I
1
0.5
0 0
Fig. 9.
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
2
A stable limit cycle around E2 where r = 1, d = 0.1, = 0.1, σ = 0, K = 2 and β = 0.91.
see that the periodic orbit expands as β increases until 0.9625039, and meets a homoclinic orbit at β = 0.9625039 (see Fig. 10, which is produced by MatCont [Dhooge et al., 2004]). Then, the homoclinic orbit is broken and all positive orbits except the infected equilibria and the stable manifolds of a positive saddle tend to the disease-free equilibrium when 0.9625039 < β < 1.2. This means that the infection dies out when β lies in this interval. Hence, higher contact coefficient can eradicate the infection. If we increase β from 1.2, two characteristic directions θ1 and θ2 of (0, 0) occur and there are bistable regions. Positive orbits in the lower part of the feasible region tend to the disease-free equilibrium, while positive orbits in the upper part
1.2 1
I
0.8 0.6
of the feasible region tend to (0, 0), which means that the infection drives the population to extinction. Furthermore, it is easily verified that the disease-free equilibrium is globally stable if 0 < β < 0.8, and the infection can be persistent for some initial values if 0.8 < β < 0.9625039. Hence, the model goes through the stages of extinction of the infection, persistence of infection, the extinction of infection, and the extinction of both infection and host population, as β increases from 0 to ∞.
3.3. Higher codimension bifurcation From the discussions above, we see that E1 and E2 coincide with E ∗ when c = 1/4. In this case, it is possible to have higher codimension bifurcations. To this end, we impose tr(J) = 0 and det(J) = 0, where J is the Jacobian matrix of (16) at an endemic equilibrium, so that we have a degenerate endemic equilibrium E ∗ . Note that y ∗ = 1/2 at the degenerate equilibrium, which implies that the determinant of the matrix J is zero. Note also that the trace of the Jacobian matrix J at the equilibrium equals zero if 1 3 −2a + H + Hσ = 0. 2 2
0.4 0.2 0
1 S
(29)
Consequently, we have E ∗ = (2a/(3 + σ), 1/2) and 0
0.5
1
1.5
S Fig. 10. Limit cycles expand as β increases from 0.91 and meet a homoclinic orbit when β = 0.9625039, where r = 1, d = 0.1, = 0.1, σ = 0, K = 2.
a1 = − a0 =
(1 + σ) (−8a + 3 + σ) , 2(3 + σ)
(1 + σ)(−8a + 3 + σ) . 4(3 + σ)
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
Theorem 17. Suppose that c = 1/4 and (29) hold.
E∗
of (16) is a cusp of codiThen the equilibrium mension 2 or higher codimension. First, we translate the equilibrium E ∗ = (2a/(3+σ), 1/2) to the origin by the transformation of X = x − 2a/(3 + σ) and Y = y − 1/2. Then we obtain dX a(1 + σ) 4a2 (1 + σ) 1 σ =− X− + X2 Y − dt 3+σ (3 + σ)2 2 2 Proof.
4aσ XY + o((X, Y )2 ), 3+σ 1 σ a (1 + σ) dY = + X+ Y + σXY dt 4 4 3+σ −
+
If Z = −(a(1 + σ)X/(3 + σ)) − 4 (a2 (1 + σ)Y/ (3 + σ)2 ), it becomes 1−σ 2 σ(3 + σ) dX =Z− X + XZ dt 2 a(1 + σ)
(30)
(3 + σ) (−3 − σ + 4a + 4aσ) 2 Z 8a2 (1 + σ)
+ o((X, Z)2 ). (3 + σ) (−3 − σ + 4a + 4aσ) XZ . P =Z+ 8a2 (1 + σ) Then (30) is transformed into 1−σ 2 dX =P− X dt 2 (3 + σ) (−4a + 4aσ + 3 + σ) XP a2 (1 + σ)
+ o((X, P )2 ), 1 + σ 2 4a(1 + σ) − 3 − σ dP = X − XP dt 8 4a + o((X, P )2 ). Set u=X+
(−3 − σ + 4a + 4aσ)X 2 , 8a
dv = dt
σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a uv 8a2 (1 + σ)
(32)
1−σ 2 u + o((u, v)2 ). 2 If σ = 1 and −
then E ∗ is a cusp of codimension 2, i.e. there is a Bogdanov–Takens bifurcation [Chow et al., 1994, Sec. 4.1]. If σ = 1 or σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a = 0,
One implication of the bifurcation of codimension 3 or higher codimension is that (16) has two limit cycles and homoclinic bifurcation under suitable perturbations.
4. Discussion
Set
+
du = v + o((u, v)2 ), dt
it follows from [Chow et al., 1994, Sec. 5.3] that E ∗ is a cusp of codimension 3 or higher codimension.
+ o((X, Z)2 ),
−
(1 + σ)X 2 . 8 Then (31) becomes v=P+
σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a = 0,
4a(1 + σ) − 3 − σ 2 Y + o((X, Y )2 ). 2(3 + σ)
1 + σ 2 4a(1 + σ) − 3 − σ dZ = X − XZ dt 8 4a
3305
(31)
Host extinction and host oscillations are reported in papers [Ebert et al., 2000; Krukonis & Schaffer 1991; Kendall et al., 1999]. Ebert et al. [2000] showed that stochastic factors are possible causes for these. The paper by Hwang and Kuang [2003] indicates that a standard incidence can lead to the extinction of a host population. This means that deterministic factors are also causes for the host extinction. In this paper, we have adopted the incidence βSI 2 /N 2 to simulate the situation where the infection probability of a susceptible individual depends linearly upon the number of exposures to infected individuals. Theorem 1 shows that the host population cannot go to extinction if β < (d + r)(1 + σ) + . This means that the population is persistent if β is small. If β > (d + r)(1 + σ) + , we have shown that the model admits saddle node bifurcations, Hopf bifurcations and codimension 2 or higher codimension
3306
W. Wang et al.
bifurcations. Numerical calculations indicate that the infection can die out even though there are two positive equilibria. The more interesting result is that the model admits bistable regions of initial states for many parameter values, in which the infection dies out if initial states lie in one region and the population and parasites coexist if initial states lie in the other region, and admits a region of initial states such that the population and parasites die out together in this region for suitable parameters. Since the extinction of the population or the infection is dependent upon initial states of the population and parasites, this seems more reasonable in reality. In contrast, either the infection dies out and the population is persistent, or both parasites and the population go to extinction together, if the standard incidence is used [Hwang and Kuang, 2003]. Hence, the correlation of infection probability of a susceptible individual to multiple exposures to infectives is a cause for bistable states. Another scenario is also interesting. With the nonlinear incidence incorporating the effect of multiple exposures, stable and unstable limit cycles can occur in model (5). Since the model with the standard incidence does not have a limit cycle (see [Hwang & Kuang, 2003]), this means that the infection probability related to multiple exposures to infected individuals is one source of oscillations of the host population. If parasites are harmless for the reproduction, i.e. σ = 1, we have shown that a Hopf bifurcation is subcritical, i.e. an unstable limit cycle is produced. For σ < 1, we find that the model has stable limit cycles. It seems that the reduction coefficient σ of birth rate of infected hosts is important for producing a stable periodic oscillation. If we decrease σ from 1 to 0, we have found in Example 13 that the model undergoes the coexistence of the population and the infection in one region of initial states, extinction of the infection, and extinction of the population and the infection in one region of initial states. Hence, medium reduction of birth rate of infected hosts seems better for infection control. The higher reduction of birth rate of infected hosts is dangerous for the population, and lower reduction increases the risk of infection. We have also found in Example 16 that the model goes through the stages of extinction of the infection, persistence of the infection in one region of initial states, extinction of the infection, and extinction of both infection and the population in one
region of initial states, as β increases from 0 to ∞. These extinction and persistence switches reveal the effects of β. High β is dangerous for the population and suitably higher β can also eradicate the infection.
References Begon, M., Bennett, M., Bowers, R. G., French, N. P., Hazel, S. M. & Turner, J. [2002] “A clarification of transmission terms in host-microparasite models: Numbers, densities and areas,” Epidemiol. Infect. 129, 147–153. Berezovskaya, F. Karev, G. & Arditi, R. [2001] “Parametric analysis of the ratio-dependent predator–prey model,” J. Math. Biol. 43, 221–246. Bouma, A., De Jong, M. C. M. & Kimman, T. G. [1995] “Transmission of pseudobabies virus within pig-populations is independent of the size of the population,” Pre Vet Med. 23, 163–172. Castillo-Chavez, C. & Yakubu, A. A. [2001] “Discretetime S-I-S models with complex dynamics,” Nonlin. Anal. 47, 4753–4762. Chow, S.-N., Li, C. Z. & Wang, D., [1994] Normal Forms and Bifurcation of Planar Vector Fields (Cambridge University Press, Cambridge). De Jong, M. C. M., Diekmann, O. & Heesterbeek, J. A. P. [1995] “How does transmission depend on population size?” in Human Infectious Diseases, Epidemic Models, ed. Mollison, D. (Cambridge University Press, Cambridge), pp. 84–94. Deredec, A. & Courchamp, F. [2003] “Extinction thresholds in host-parasite dynamics,” Ann. Zool. Fennici 40, 115–130. Derrick, W. R. & van den Driessche, P. [1993] “An infection transmission model in a nonconstant population,” J. Math. Biol. 31, 495–512. Derrick, W. R. & van den Driessche, P. [2003] “Homoclinic orbits in an infection transmission model,” Discr. Contin. Dyn. Syst. Ser. B. 3, 299–309. Dhooge, A., Govaerts, W. & Kuznetsov, Y. A. [2004] Limit cycles and their bifurcations in MatCont. http://allserv.UGent.be/ajdhooge. Diekmann, O. & Heesterbeek, J. A. P. [2000] Mathematical Epidemiology of Infectious Diseases. Model Building, Analysis and Interpretation, Wiley Series in Mathematical and Computational Biology (John Wiley & Sons, Ltd., Chichester). Ebert, D., Lipsitch, M. & Mangin, K. L. [2000] “The effect of parasites on host population density and extinction: Experimental epidemiology with Daphnia and six microparasites,” Amer. Naturalist 156, 459– 477. Gao, L. & Hethcote, H. W. [1992] “Disease transmission models with density-dependent demographics,” J. Math. Biol. 30, 717–731.
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
Grenfell, B. T. & Dobson, A. P. (eds.) [1995] Ecology of Disease in Natural Populations (Cambridge University Press, Cambridge). Guckenheimer, J. & Holmes, P. J. [1996] Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer-Verlag, NY). Hethcote, H. W. [2000] “The mathematics of infectious diseases,” SIAM Rev. 42, 599–653. Hwang, T. W. & Kuang, Y. [2003] “Deterministic extinction effect of parasites on host populations,” J. Math. Biol. 46, 17–30. Kendall, B. E., Briggs, C. J., Murdoch, W. W., Turchin, P., Ellner, S. P., McCauley, E., Nisbet, R. M. & Wood, S. N. [1999] “Why do populations cycle? A synthesis of statistical and mechanistic modelling approaches,” Ecology 80, 1789–1805. Krukonis, G. & Schaffer, W. M. [1991] “Population cycles in mammals and birds: Does periodicity scale with body size?” J. Theor. Biol. 148, 469–493. Liu, W. M., Levin, S. A. & Iwasa, Y. [1986] “Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models,” J. Math. Biol. 23, 187–204. Liu, W. M., Hethcote, H. W. & Levin, S. A. [1987] “Dynamical behavior of epidemiological models with nonlinear incidence rates,” J. Math. Biol. 25, 359–380. Lizana, M. & Rivero, J. [1996] “Multiparametric bifurcations for a model in epidemiology,” J. Math. Biol. 35, 21–36. May, R. M. & Anderson, R. M. [1979] “Population biology of infectious diseases. II,” Nature 280, 455–461.
3307
McCallum, H. [1995] “Modelling wildlife-parasite interactions to help plan and interpret field studies,” Wildlife Res. 22, 21–29. Nold, A. [1980] “Heterogeneity in disease-transmission modeling,” Math. Biosci. 52, 227–240. Perko, L. [1996] Differential Equations and Dynamical Systems (Springer-Verlag, NY). Polking, J. C., http://math.rice.edu/dfield/6.5. Ruan, S. & Wang, W. [2003] “Dynamical behavior of an epidemic model with a nonlinear incidence rate,” J. Diff. Eqs. 188, 135–163. van den Driessche, P. & Watmough, J. [2002] “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Math. Biosci. 180, 29–48. van Nouhuys, S. & Ehrnsten, J. [2004] “Wasp behavior that leads to uniform parasitism of a host available only a few hours per year,” Behav. Ecol., in press. Wang, W. & Ruan, S. [2004] “Simulating the SARS outbreak in Beijing with limited data,” J. Theor. Biol. 227, 369–379. Wang, W. & Zhao, X.-Q. [2004] “An epidemic model in a patchy environment,” Math. Biosci. 190, 97– 112. Xiao, D. & Ruan, S. [2001] “Global dynamics of a ratiodependent predator–prey system,” J. Math. Biol. 43, 268–290. Zhang, Z., Ding, T., Huang, W. & Dong, Z. [1991] Qualitative Theory of Differential Equations, Trans Mathematical Monographs, Vol. 101 (Amer. Math. Soc., Providence).
BIFURCATIONS IN A HOST-PARASITE MODEL WITH NONLINEAR INCIDENCE WENDI WANG∗ Department of Mathematics, Southwest Normal University, Chongqing 400715, P. R. China YI LI and H. W. HETHCOTE Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA Received May 20, 2005; Revised October 26, 2005 A host-parasite model is proposed that incorporates a nonlinear incidence rate. Under the influence of multiple infectious attacks, the model admits bistable regions such that the infection dies out if initial states lie in one region, and the population and parasites coexist if initial states lie in the other region. It is also found that parasites can drive the population to extinction for suitable parameters. It is verified that the model has a saddle-node bifurcation, Hopf bifurcations and a cusp of codimension 2 or higher codimension. Stable limit cycles and unstable limit cycles are examined as the infection-reduced reproduction rate varies. It is shown that the model goes through the change of stages of infection extinction, infection persistence, infection extinction, and the extinction of both parasites and the population as the contact coefficient increases. Keywords: Nonlinear incidence; host extinction; periodic solutions; homoclinic orbits; bistable.
1. Introduction Classical epidemic models are extended in many ways to understand mechanisms of disease transmissions (see, e.g., [Hethcote, 2000; Wang & Zhao, 2004; Wang & Ruan, 2004]). An interesting topic in the study of infectious diseases is to understand how parasites regulate host populations (see [May & Anderson, 1979; De Jong et al., 1995; Grenfell et al., 1995; McCallum, 1995]). Recently, with the aim of understanding how six microparasites regulate Daphnia populations and drive the populations to extinction, Ebert et al. [2000] formulated the following epidemiological microparasite ∗
model with horizontal transmission: S+I dS = r(S + σI) 1 − − dS − λSI , dt K
(1)
dI = λIS − ( + d)I, dt where S(t) and I(t) represent the densities of uninfected (susceptible) and infected (infective) hosts at time t, respectively, r is the per capita growth rate of uninfected hosts, σ is the relative fecundity of an infected host with 0 ≤ σ ≤ 1, K is the carrying capacity of the environment for the host population, d is the parasite mortality rate, λ is the infection
Research is supported by the National Natural Science Fund of China (No. 10571143). 3291
3292
W. Wang et al.
rate coefficient, and is the infection-induced death rate. Model (1) uses the classical mass action incidence λIS and has the usual asymptotic behavior. If there is no endemic equilibrium, the infection free equilibrium is globally stable. If it does have an infected equilibrium, this steady state is globally stable. Thus this model fails to explain the observed rich outcomes that depend on parameter values and initial population levels. To explain why the host population can be driven to extinction, a carefully designed stochastic simulation of model (1) was conducted in [Ebert et al., 2000]. The simulations indicated that extinction of the host is possible in some parameter regions. This means that stochastic factors are possible causes of the extinction of the host population. The mass-action principle stems from the theory of the kinetics of chemical reactions. It implies that individuals, like molecules, would be mixed homogeneously in space and every encounter would have the same occurrence probability. This may be applicable at low host population densities. For large densities, the number of new encounters made by a single individual per unit time cannot increase linearly with population sizes, either because overcrowding reduces the movements of individuals, or because a related increase of occupied habitat prevents its total occupation by individuals. In this case, a standard incidence λSI/(S + I) would be a suitable option (see [Hethcote, 2000] and [Deredec et al., 2003]). Diekmann and Heesterbeek [2000] gave a more general incidence C(N )SI/N with N = S +I, where C(N ) is the number of individuals that are engaged in contacts at any time. This incidence is suitable for any population size, and reduced to the mass action law if C(N ) is linearly dependent on N at low population densities, or reduced to the standard incidence if C(N ) takes a constant value at high population densities. However, the standard incidence may also apply to low population densities due to adaptive behavior of populations. It is observed that larger populations may occupy larger areas [Begon et al., 2002]. This phenomenon seems more common for wild populations. Based upon these observations, [De Jong et al., 1995] and [Bouma et al., 1995] suggested that each individual occupies a characteristic area such that the area occupied by a population is directly proportional to the size of that population. In this case, the incidence at low host population size may also be represented or approximated by the standard incidence.
Hwang and Kuang [2003] changed the mass action incidence in (1) to the standard incidence λSI/(S + I). Significantly, this revised parasite-host model can exhibit the observed parasite-induced host extinction. This means that the extinction of the host population may be due to deterministic factors, instead of stochastic factors. Since the mass action incidence implies that individuals of the population occupy a fixed area and the standard incidence implies that the occupied area is proportional to its population size, alternately, the extinction of the host population may be from the behavior changes of the population from the movement in a fixed area to the movement in an area that is proportional to its population size. However, periodic fluctuations of the population and outcomes that are dependent on initial population levels are not predicted in the model with the standard incidence, as pointed out in their paper [Hwang & Kuang, 2003]. Note that the standard incidence can be further written as β0 ωSI/(S + I) where ω is the encounter rate and β0 is the probability that a susceptible individual becomes infected in an encounter by an infectious individual. The standard incidence implies that the infection probability is invariant for all sizes of infectious individuals and susceptibility of the individuals. However, there is evidence that infections are limited by host physiological resistance and behavioral defense of the immune system [Nouhuys et al., 2004]. In these cases, a susceptible is more vulnerable under multiple attacks such that multiple attacks by infectious members may increase the infection probability of the susceptible. Thus, it is reasonable to assume that β0 is an increasing function of the number of infectious attacks ωI/(S + I) per unit time. To explore the influence of this infection probability as a function of multiple attacks on the dynamics of hostparasite evolutions, in this paper we consider only the simplest case that the infection probability is proportional to the number of infectious attacks ωI/(S + I). Then we have β0 = β1 ωI/(S + I), where β1 is a proportional constant. As a consequence, the incidence becomes βSI 2/(S + I)2 with β = β1 ω 2 . Indeed, this incidence is a special case of the infection force βG(I/N ) of proportional mixing which was proposed and used by Nold [1980], Derrick and van den Driessche [1993], Castillo-Chavez and Yakubu [2001] and Derrick and van den Driessche [2003]. If a population size is constant, this type of nonlinear incidences was studied by [Liu et al.,
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
1986, Liu et al., 1987; Lizana & Rivero, 1996; Ruan & Wang, 2003]. Based upon the incidence and the modeling approach for the endemic models with the logistic demographic structure in [Gao & Hethcote, 1992], we formulate the following SI endemic model: χrN dN = b− [S + σI] dt K (1 − χ)rN N − I − d+ K N N − ( + b(1 − σ))I =r 1− K χrIN K χrN = b− [S + σI] K I2 (1 − χ)rN S − βS 2 , − d+ K N I2 (1 − χ)rN = βS 2 − + d + I, N K + (1 − σ)
dS dt
dI dt
(2)
where N = S + I is the total host population density, b − χrN/K is the per capita birth rate for susceptibles, and σ is the infection-induced reduction in the birth rate, so that [b − χrN/K]σ is the per capita birth rate for infectives. Here d+(1−χ)rN/K is the natural per capita death rate, r = b − d > 0 is the per capita net growth rate, K is the environmental carrying capacity, χ is the convex combination constant with 0 ≤ χ ≤ 1, and is the per capita infection-related death rate. Note that the population dynamics are governed by the logistic differential equation, when there is no infection (I = 0). For 0 < χ < 1, the per capita birth rate b − χrN/K is a decreasing function of N and the natural per capita death rate d + (1 − χ)rN/K is an increasing function of N , which is consistent with densitydependent population dynamics. When χ = 1, the density dependence affects the per capita birth rate, but does not affect the per capita death rate. In order to get an SI epidemic model comparable to those considered by Elbert et al. [2000], Hwang and Kuang [2003], we set χ = 1, so that the model (2) becomes rN I2 dS = b− (S + σI) − dS − βS 2 , dt K N (3) 2 I dI = βS 2 − ( + d)I. dt N
3293
This formulation of an SI endemic model with the density-dependent growth is slightly different from that in model (1), in which dN N = r[N − (1 − σ)I] 1 − − dN − I. dt K Thus, in the absence of the disease, model (1) implies rN dN = r−d− N dt K so that N → (r − d)K/r as t → ∞. The advantage of the model (3) is that the original meanings of parameters r and K in the logistic differential equation are maintained, in the sense that r is the intrinsic growth rate of the population and N → K as t → ∞ in the absence of the disease. For biological reasons, we will restrict out attention on the region {(S, I) : 0 ≤ S + I ≤ K}. Note that (3) has the singularity at (0, 0). We make a time scale change dt = N 2 dτ such that (3) is equivalent to the following system in the interior of the first quadrant R2 : dS = (r + 2bσ − β)SI 2 + (2r + bσ)S 2 I + rS 3 dt + bσI 3 −
rS 4 rσI 4 r(3 + σ)S 3 I − − K K K
r(1 + σ)S 2I 2 r(3σ + 1)SI 3 − , −3 K K
(4)
dI = I[βSI − ( + d)(S + I)2 ], dt where t represents τ for the convenience of notations. To be concise in analysis, let us scale (4) by x = rS/(βK), y = rI/(βK) and ξ = β 3 K 2 t/r 2 . If t is used to represent ξ, we obtain dx = ax3 + (2a + σh)x2 y + (a + 2σh − 1)xy 2 dt + σhy 3 − x4 − (3 + σ)x3y − 3(1 + σ)x2 y 2 − (3σ + 1)xy 3 − σy 4 ,
(5)
dy = y(xy − c(x + y)2 ), dt where a=
r , β
h=
d+r , β
c=
d+ . β
In this paper, we assume that r, β, d and are positive constants, and σ is a constant in [0, 1]. Hence, we
3294
W. Wang et al.
have h > a and a + c − hσ > 0. The feasible region of (5) is {(x, y) : x ≥ 0, y ≥ 0, 0 ≤ x + y ≤ a}. The organization of this paper is as follows. In the next section, we analyze dynamical behavior of (5) around (0, 0) and give conditions under which the infection cannot drive the population to extinction, and conditions under which the infection drives the population to extinction. In Sec. 3, we present the analysis of bifurcation structures including Hopf bifurcations and higher codimension bifurcations, and show the existence of limit cycles. We conclude with a discussion of the results in Sec. 4.
2. Extinctions In this section, we analyze the existence and stability of boundary equilibria of (5) and present conditions under which the infection drives the population to extinction and conditions under which the infection cannot drive the population to extinction. First, (a, 0) is an infection free equilibrium. It is always asymptotically stable. This means that the outcome of infection is unsuccessful if the invasion of parasites is small. To be endemic, the intensity of invasion must be higher than a threshold. Secondly, E0 = (0, 0) is a trivial equilibrium. This corresponds to the extinction of both infected hosts and the host population. Since E0 is nonhyperbolic, we introduce the polar coordinates x = R cos θ, y = R sin θ to obtain dR = R3 (V (θ) + o(1)), dt
(6)
dθ = −R2 sin θ(cos θ + sin θ)(U (θ) + o(1)), dt where U (θ) = (a + c)cos2 θ + (a + c + hσ − 1)cos θ sin θ + hσ sin2 θ, V (θ) = a(cos(θ))2 (sin(θ) + cos(θ))2 + hσ sin(θ)cos(θ)(sin(θ) + cos(θ))2 − c(sin(θ))2 (sin(θ) + cos(θ))2 + (sin(θ))2 cos(θ)(− cos(θ) + sin(θ)). Note that the first quadrant R2 is positively invariant. By the qualitative theory of differential equations [Zhang et al., 1991] (see also [Xiao & Ruan, 2001; Berezovskaya et al., 2001] for population models), an orbit of (5), if it tends to the origin, must tend to it along a characteristic direction. θ = 0 is a trivial characteristic direction, where infected hosts die out and the uninfected
host population tends to its carrying capacity K. We now look for other characteristic directions in (0, π/2). If z = tan θ, by the definition of U we see that z satisfies F1 (z) := hσz 2 + (a + c + hσ − 1)z + a + c = 0. (7) Thus, there is no characteristic direction in (0, π/2) if a + c + hσ − 1 > 0 or
a + c + hσ − 1 < 0, (a + c + hσ − 1)2 < 4hσ(a + c).
(8)
(9)
If (8) or (9) is satisfied, since V (0) = a > 0 and U (0) = a + c, by Theorems 3.5 and 3.7 of [Zhang et al., 1991], we have the following theorem: Theorem 1. If (8) or (9) is satisfied, then there is
an 1 > 0 and an R1 > 0 such that (5) admits a unique orbit in {(θ, R) : 0 ≤ θ < 1 , 0 < R < R1 } that tends to (0, 0) along θ = 0 as t → −∞. Biological implication of Theorem 1 is that the population and parasites cannot die out together because (0, 0) is a repeller. Note that (7) has two positive solutions z1 and z2 with z1 < z2 if a + c + hσ − 1 < 0, (10) (a + c + hσ − 1)2 > 4hσ(a + c). This means that we have two other characteristic directions 0 < θ1 < θ2 except for θ = 0. In order to find dynamical behaviors of (5) along these two directions, we apply the Briot–Bouquet transformation: x = x,
y = ux,
dτ = x2 dt
to transform (5) into dx = x[a − x + (2a + hσ)u − (3 + σ)xu dt + (a + 2hσ − 1)u2 − 3(1 + σ)xu2 + hσu3 − (3σ + 1)xu3 − σxu4 ], du = −u(u + 1)[a + c − x dt
(11)
+ (a + c + hσ − 1)u − (2 + σ)xu + hσu2 − (1 + 2σ)xu2 − σxu3 ], where t is used to represent τ . We consider (11) in the region x ≥ 0, u ≥ 0 due to the biological
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
background. Under the condition (10), we see that (11) admits equilibria (0, 0), (0, z1 ) and (0, z2 ). It is easy to see that (0, 0) is a saddle. To analyze equilibrium (0, z1 ), we make the change of variables x1 = x, x2 = u − z1 to transform (11) to dx1 = [a + (2a + hσ)z1 + (−1 + a + 2hσ)z12 dt + hσz13 ]x1 + o((x1 , x2 )), dx2 = z1 (1 + z1 )3 (σz1 + 1)x1 + [−a − c dt + (−4a − 2hσ − 4c + 2)z1 (−6hσ + 3 − 3a − 3c)z12 − 4hσz13 ]x2 + o((x1 , x2 )),
(12)
Using (7), (12) can be simplified to dx1 = g1 (z1 )x1 + o((x1 , x2 )), dt dx2 = z1 (1 + z1 )3 (σz1 + 1)x1 dt + g2 (z1 )x2 + o((x1 , x2 )),
Further, (0, z1 ) of (11) is a saddle if hσ < h3 σ < c(1 + (a/(c − 1))) and F2 (h) > 0. Note that A0 := a + c − hσ > 0. Condition (i) means
Proof.
A1 := −cσh − c + ca + c2 + hσ ≥ 0. Hence, g1 (z1 ) = A1 z1 + A0 > 0. If (ii) holds, then A1 < 0. Since h ≥ h3 , we have F0 (h) < 0. By direct calculations, we see that F0 (h) has the same sign as A0 a + c + hσ − 1 . − A1 2hσ Hence, A0 > 0. g1 (z1 ) = A1 z1 − −A1 When (iii) holds, we have a + c + hσ − 1 A0 . > A1 2hσ
(13)
where g1 (z1 ) = c(a − hσ + c) + (−cσh − c + ca + c2 + hσ)z1 , g2 (z1 ) = −(a + c)(a − hσ + c − 1) + (hσ(c + a + 1) − (a + c − 1)2 )z1 . The stability of (0, z1 ) is determined by the signs of gi (z1 ). Let us consider (a + c + hσ − 1)2 − 4 hσ(a + c) = 0.
3295
(14)
This equation in h has two solutions 0 < h1 < h2 . It is easy to verify that (10) implies 0 < h < h1 . Next, we consider F0 (h) := h2 σ 2 (c + 1) − hσ(−a + 2c2 − c + 2ac + 1) + c(a + c − 1)2 = 0. (15) This equation in h has two solutions 0 < h3 < h4 with h1 < h4 < h2 . Set F2 (h) = ch2 σ 2 − hσ(2c2 − a + 2ac) + c(c2 + 2ac + a2 − a).
Note that A0 hσ = F1 −A1 (−σhc − c + ac + c2 + hσ)2 × F2 (h) < 0. Thus, A0 /(−A1 ) > z1 . Consequently, we have g1 (z1 ) > 0. Now, we consider the sign of g2 (z1 ). Since a + c + hσ < 1, we have B0 := −(a + c)(a − hσ + c − 1) > 0. Set B1 := hσ(c + a + 1) − (a + c − 1)2 . Using the fact that h < h1 , we obtain hσ < (a + c − 1)2 /(a + c + 1), which means B1 < 0. Then by similar arguments to above, we obtain a + c + hσ − 1 B0 , F2 (h) < 0, > B1 2hσ where h < h1 is used. Thus, B0 /(−B1 ) > z1 . It follows that B0 > 0. g2 (z1 ) = B1 z1 − −B1
Theorem 2. Let (10) hold. Then (0, z1 ) of (11) is an unstable node if one of the following conditions is satisfied:
In summary, we have verified that g1 (z1 ) > 0 and g2 (z1 ) > 0 if any of (i)–(iii) is satisfied. Thus, (0, z1 ) is an unstable node in this case. The conclusion that (0, z1 ) is a saddle can be proved in a similar way.
(i) c(1 + (a/(c − 1))) < hσ; (ii) h3 σ ≤ hσ < c(1 + (a/(c − 1))); (iii) hσ < h3 σ < c(1 + (a/(c − 1))) and F2 (h) < 0.
By similar discussions to those in the proof of Theorem 2, we can obtain the following results for the stability of (0, z2 ).
3296
W. Wang et al.
0.4
0.35
0.3
I
0.25
0.2
0.15
0.1
0.05
0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
S
0.1
0.09
0.08
0.07
I
0.06
0.05
0.04
0.03
0.02
0.01
0 0
0.05
0.1
0.15
S
Fig. 1.
The figures of two typical ways of extinction of both host population and parasites.
Theorem 3. Let (10) hold. Then (0, z2 ) of (11) is a saddle if one of the following conditions is satisfied:
(i) c(1 + (a/(c − 1))) < hσ; (ii) h3 σ < hσ < c(1 + (a/(c − 1))) and F2 (h) > 0. Further, (0, z2 ) of (11) is a stable node if one of the following holds: (1) h3 σ < hσ < c(1 + (a/(c − 1))) and F2 (h) < 0; (2) hσ < h3 σ < c(1 + (a/(c − 1))). From Theorems 2 and 3, there are three types of dynamical behaviors around (0, 0) for (5) when the characteristic directions θ1 and θ2 exist. First, if (0, z1 ) is an unstable node and (0, z2 ) is a saddle, then (0, 0) of (5) is a repeller so that the host population and parasites cannot die out together.
Secondly, if (0, z1 ) is a saddle and (0, z2 ) is a stable node, the region in the first quadrant around (0, 0) is split into two parts, one of which is an attracting parabolic sector and the other is a hyperbolic sector. The host population and parasites go to extinction in the attracting parabolic sector. Lastly, if (0, z1 ) is an unstable node and (0, z2 ) is a stable node, the region in the first quadrant around (0, 0) is split into two parts, one of which is a repelling parabolic sector and the other is an elliptic sector. The host population and parasites go to extinction in the elliptic sector (see Fig. 1).
3. Bifurcations 3.1. Positive equilibria We now consider positive equilibria of (5). Roughly speaking, the existence of these equilibria is
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
necessary for the coexistence of the population and parasites. First, we perform a topological transformation to (5) so that positive equilibria can be easily found. Set u = x + y, v = y/(x + y) and dτ = x2 dt. Substituting them into (5) and then replacing u by x, v by y, τ by t for the simplicity of notations, we obtain dx = x(a − x + (hσ − c − a)y + (1 − σ)xy), dt dy = y(−c − a + x + (1 + c + a − hσ)y dt − (1 − σ)xy − y 2 ).
(16)
Evidently, we should consider this system in the region D = {(x, y) : 0 < x < a, 0 < y < 1}. The isoclines of (16) in D are given by x= x=
a − (a − hσ + c)y = f1 (y), 1 − (1 − σ) y
a + c − (1 + c + a − hσ)y + 1 − (1 − σ)y
y2
= f2 (y).
(17) (18)
(19)
If hσ < c, we have a/(a + c − hσ) < 1. Thus, it suffices to consider solutions of (19) in (0, a/(a + c − hσ)]. Note that f3 (0) = c > 0 and F2 (h) a = . f3 a + c − hσ (a + c − hσ)2
x(−σy + y − 1) J := y(σy − y + 1)
If c ≥ hσ + a, then a/(a + c − hσ) ≤ 1/2. Note that f3 is a decreasing function on (0, a/(a + c − hσ)). Consequently, (19) has a unique solution on (0, a/(a + c − hσ)) if F2 (h) < 0, and has no solution on (0, a/(a + c − hσ)) if F2 (h) > 0. Similarly, if hσ < c < hσ + a, then (19) has a unique solution on (0, a/(a + c − hσ)) if c < 1/4 and F2 (h) < 0, has two solutions on (0, a/(a + c − hσ)) if c < 1/4 and F2 (h) > 0, and has no solution on (0, a/(a+c−hσ)) if c > 1/4. For the case where c ≤ hσ, since f3 (1) = c > 0, it follows that (19) has two positive solutions on (0, 1) if c < 1/4, has a unique positive solution on (0, 1) if c = 1/4 and has no solution on (0, 1) if c > 1/4. In summary, we have the following results for positive equilibria of (16). Theorem 4. (16) has a unique positive equilibrium
E ∗ if one of the following conditions holds:
To obtain a positive equilibrium in D, we should have y < a/(a + c − hσ). By (17) and (18), we obtain f3 (y) := y 2 − y + c = 0.
3297
(i) c ≥ hσ + a and F2 (h) < 0; (ii) hσ < c < min{hσ + a, 1/4} and F2 (h) < 0; (iii) c = 1/4 < hσ. Further, (16) has two positive equilibria E1 = (x1 , y1 ) and E2 = (x2 , y2 ) where y1 < 1/2, y2 > 1/2 and xi = f1 (yi ), i = 1, 2 if one of the following conditions is satisfied: (1) hσ ≤ c < min{hσ + a, 1/4} and F2 (h) > 0; (2) c < min{hσ, 1/4}. Now, we consider the stability of positive equilibria. The Jacobian matrix at a positive equilibrium (x, y) is
−x(xσ − x + c + a − hσ) . y(1 − 2y + xσ − x + c + a − hσ)
Its characteristic equation is λ2 + (σyh − ay − yc + x − y + 2y 2 )λ + yx(2y − 1)(σy − y + 1) = 0. Thus, det(J) < 0 if y < 1/2 and det(J) > 0 if y > 1/2. Consequently, we can state. Theorem 5. If (16) has a unique endemic equilib-
rium E ∗ and (iii) in Theorem 4 does not hold, this endemic equilibrium is a saddle. Further, if (16) has two endemic equilibria, then E1 is a saddle, E2 is a node or focus. One implication of Theorems 4 and 5 is the existence of a saddle node bifurcation. For example, if
we fix d = 0.2, = 0.1, r = 0.16, σ = 0.6, then a = 0.16/β, h = 0.36/β, c = 0.3/β. In this case, there are two positive equilibria when 1.2 < β < 1.3289, one of which is a saddle and the other one is a node or a focus (see Fig. 2). To indicate another way of occurrence of positive equilibria, we set d = 0.2, = 0.1, r = 0.16, β = 1.22. Then a = 0.1311, h = 0.2951, c = 0.2459. It is easy to see that (16) does not have a positive equilibrium if 0 ≤ σ < 0.25, has a unique positive equilibrium if 0.25 < σ < 0.4898. When σ increases from 0.4898, (0, 0) becomes a repeller and the second positive equilibrium emerges.
W. Wang et al.
3298
If a1 < 0 and a0 ≤ 0, then the trace is negative. Therefore, E2 is stable. If a1 < 0 and a0 > 0, we have a0 tr 1 = a1 y2 − . (20) −a1
0.07 0.06 0.05 0.04
x
Evidently, tr 1 < 0 if 0.03
2a0 + a1 < 0. 0.02 0.01 0 −0.01 1.18
1.2
Fig. 2.
1.22
1.24
1.26 β
1.28
1.3
1.32
1.34
The graph of a saddle node bifurcation.
We now determine the stability of E2 = (x2 , y2 ). Set H = a + c − hσ and a1 = H(1 + σ) − 2c − σ + 2cσ a0 = H(−1 + c − cσ) − hσ + 2c + cσ, r1 = (c2 (1 − σ)2 + σ(4c − 1))H 2 + H(hσ(−1 + 2c)(−1 + σ) + σ − 4cσ − 2c2 + 2c2σ) + h2σ 2 − hσ(2c + 4cσ − σ) + c(−2σ − cσ 2 + 10cσ + 4c2 (1 − σ)2 ). Theorem 6. Assume that (16) has two endemic equilibria. Then E2 is asymptotically stable if any one of the following conditions is satisfied
(1) a1 < 0 and 2a0 + a1 < 0; (2) a1 < 0, 2a0 + a1 > 0, r1 < 0; (3) a1 > 0, 2a0 + a1 < 0 and r1 > 0. Further, E2 is unstable if any one of the following conditions is satisfied (i) a1 > 0 and 2a0 + a1 > 0; (ii) a1 > 0, 2a0 + a1 < 0 and r1 < 0; (iii) a1 < 0, 2a0 + a1 > 0 and r1 > 0. Note that the determinant of the Jacobian matrix J is positive at E2 . It suffices to consider the sign of the trace of J. Note that the sign of the trace is determined by Proof.
(−x + y − σhy2 + ay2 +
cy2 − 2y22 ).
Using (19) and (17), we see that its sign is determined by tr 1 := a1 y2 + a0 .
(21)
Hence, E2 is stable if condition (1) holds. By direct calculations, we see that the sign of f3 (−a0 /a1 ) is determined by r1 . If a1 < 0, a0 > 0 and 2a0 + a1 > 0, which implies −a0 /a1 > 1/2, it follows from r1 < 0 that −a0 /a1 < y2 , which leads to tr 1 < 0. Hence, E2 is stable if condition (2) holds. Similarly, E2 is stable if condition (3) is satisfied. The proofs that E2 is unstable when one of the conditions (i)–(iii) holds are similar to the above discussions. At this stage, we are able to indicate certain global behaviors of (5). First, we consider the case where (5) does not have a positive equilibrium. Then there are two types of dynamical behaviors. If (0, 0) is a repeller, the infection-free equilibrium is globally stable, which is stated in Theorems 7 and 8. If there is an attracting sector, two split regions occur where the host population and parasites die out in one region and the infection dies out in the other region. This case is described in Theorem 9. Theorem 7. Let (8) or (9) hold. Then the infection-
free equilibrium (a, 0) is globally stable if one of the following conditions is satisfied: (i) c > hσ + a and F2 (h) > 0 holds; (ii) c < hσ + a and c > 1/4. Theorem 1 implies that (0, 0) repels positive orbits if (8) or (9) holds. By the arguments preceding Theorem 4, we see that (5) does not have a positive equilibrium if (i) or (ii) holds. It is easy to verify that positive solutions of (5) are ultimately bounded. It follows that positive solutions of (5) tend to the infection-free equilibrium. Since it is always stable, the global stability follows.
Proof.
Similarly, we can obtain the following results: Theorem 8. Let (10) hold. Assume that one of the
conditions (i)–(iii) in Theorem 2 is satisfied and that one of the conditions (i) and (ii) in Theorem 3 is satisfied. Then the infection-free equilibrium (a, 0) is globally stable if one of the following conditions
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
is satisfied:
all positive orbits of (5) initiating from the other region tend to (0, 0) along the characteristic direction θ2 as t tends to infinity.
(i) c > hσ + a and F2 (h) > 0 holds; (ii) c < hσ + a and c > 1/4. Theorem 9. Let (10) hold. Assume further that
one of the conditions (1) and (2) in Theorem 3 is satisfied. If one of the following conditions:
Example 11. Let us consider the differential equa-
tions in (3). r(S + I) dS = (S + σI) b − − dS dt K 2 I S, −β S+I 2 I dI =β S − ( + d)I. dt S+I
(i) c > hσ + a and F2 (h) > 0 holds, (ii) c < hσ + a and c > 1/4, is satisfied, then there is a positive orbit in (5) which tends to (0, 0) along the characteristic direction θ1 as time tends to ∞ or −∞, and which splits the interior of (5) into two regions such that all positive orbits of (5) initiating from one region tend to (a, 0) as t tend to infinity, and all positive orbits of (5) initiating from the other region tend to (0, 0) along the characteristic direction θ2 as t tends to infinity. We have similar results if there is a unique positive equilibrium which is nondegenerate. Theorem 10. Let (10) hold. Assume further that
one of the conditions (1) and (2) in Theorem 3 is satisfied. If one of the following conditions: (i) c ≥ hσ + a and F2 (h) < 0; (ii) hσ < c < min{hσ + a, 1/4} and F2 (h) < 0; is satisfied, then the stable manifolds of E ∗ splits the interior of feasible region of (5) into two regions such that all positive orbits of (5) initiating from one region tend to (a, 0) as t tend to infinity, and
0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 3.
0.2
0.4
0.6
0.8
(22)
We fix r = 0.5, d = 0.4, σ = 0.1, = 0.5. After we transform (22) into (5), we obtain a = 0.5/β, h = 0.9/β and c = 0.9/β. When 0 < β < 2.1999, Theorem 7 is satisfied. Thus, the infectionfree equilibrium is globally stable. When 2.1999 < β < 3.8136, Theorem 9 is satisfied where there is no endemic equilibrium. When β > 3.8136, Theorem 10 is satisfied where there is a unique endemic equilibrium. Thus, in the last two cases, there are bistable regions and the outcome of evolution of the population and parasites depends on their initial conditions (see Figs. 3 and 4. This figure and others in the following are produced by the package PPLANE6 [Polking, 2003] unless stated otherwise). The infection dies out and the population tends to its carrying capacity in the lower region, while both the population and parasites die out in the upper region.
1
0
3299
1 S
1.2
1.4
1.6
1.8
2
The graph of bistable regions where r = 0.5, d = 0.4, σ = 0.1, = 0.5, K = 2 and β = 2.5.
3300
W. Wang et al. 1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 4. The stable manifolds of E ∗ split the feasible region into two parts. The lower part is a basin of the infection-free equilibrium and the upper part is a basin of (0, 0), where r = 0.5, d = 0.4, σ = 0.1, = 0.5, K = 2 and β = 4.5.
We should address that there is an elliptic sector in the case where Theorem 10 is satisfied (see Fig. 4). When two positive equilibria occur, global structure of (5) is more complicated because periodic solutions and homoclinic orbits may occur. To obtain information of periodic solutions of (16), we consider its Hopf bifurcations.
Theorem 12. Let condition (1) or condition (2) in Theorem 4 hold and r1 = 0. Assume that a1 < 0, 2a0 + a1 > 0 or a1 > 0, 2a0 + a1 < 0. Then there is a stable periodic orbit in (16) near the Hopf bifurcation point when ξ > 0, and there is an unstable periodic orbit in (16) near the Hopf bifurcation point when ξ < 0.
In order to determine the direction of a Hopf bifurcation, let us make the transformation of X = x − x2 , Y = y − y2 . Then (16) becomes
Proof.
3.2. Hopf bifurcation A Hopf bifurcation may occur when c < min{hσ + a, 1/4} and r1 = 0. In order to know the direction of a Hopf bifurcation, we set ξ := ((a1 + 4a0 )(1 − σ) + 2a1 )a1 H − + (a21 + 4a1 a0 + 6a20 )(1 − σ).
a21
dX = a11 X + a12 Y + g1 (X, Y ), dt
(23)
dY = a21 X + a22 Y + g2 (X, Y ), dt
where a11 = a12 =
a0 (a0 (1 − σ) + a1 )(Ha1 + a1 + 2a0 ) , a31
(Ha 21 + (1 − σ)(a1 a0 + a1 a0 H + 2a20 ))(Ha 1 + a1 + 2a0 )a0 , a41
a21 = −
a0 (a0 (1 − σ) + a1 ) , a21
a22 = −a11 and gi , i = 1, 2, are terms of higher order. Now, let Z = a11 X + a12 Y to obtain Z − a11X dZ Z − a11X Z − a11X dX = Z + g1 X, = −a2 X + a11 g1 X, , + a12 g2 X, , dt a12 dt a12 a12 where a2 = −
(a1 + 2a0 )(a0 (1 − σ) + a1 )(Ha1 + a1 + 2a0 )a20 . a51
(24)
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
3301
√ Set u = −X, v = Z/ a2 . Then (24) becomes √ du = − a2 v + G1 (u, v), dt where
Set
dv √ = a2 u + G2 (u, v), dt
(25)
√ v a2 + a11 u , G1 (u, v) = −g1 −u, a12 √ √ v a2 + a11 u v a2 + a11 u a11 g1 −u, + a12 g2 −u, a12 a12 . G2 (u, v) = √ a2 2 2 ∂ 3 G2 ∂ 3 G1 ∂ 3 G2 1 ∂ G1 ∂ G1 ∂ 2 G1 1 ∂ 3 G1 + + 2 + + + √ Φ(u, v) = 16 ∂u3 ∂u∂v 2 ∂u ∂v ∂v 3 16 a2 ∂u∂v ∂u2 ∂v 2 ∂ 2 G1 ∂ 2 G2 ∂ 2 G1 ∂ 2 G2 ∂ 2 G2 ∂ 2 G2 ∂ 2 G2 + + − . − ∂u∂v ∂u2 ∂v 2 ∂u2 ∂u2 ∂v 2 ∂v 2
By the results in [Guckenheimer & Holmes, 1996] or [Perko, 1996], the direction of the Hopf bifurcation is determined by the sign of Φ(0, 0). The expression of Φ(0, 0) is very complicated. However, using the facts that a11 + a22 = 0 and y22 − y2 + c = 0, we have c=−
a0 (a1 + a0 ) , a21
a=−
a0 ((a1 a0 H + a1 a0 + 2a20 )(1 − σ) + 2Ha21 + 2a0 a1 + a21 ) . a31
As a consequence, we see that the sign of Φ(0, 0) is determined by ξ where ξ : = ((a1 + 4a0 )(1 − σ) + 2a1 )a1H − a21 + (a21 + 4a1 a0 + 6a20 )(1 − σ). The conclusion of Theorem 12 now follows from [Guckenheimer & Holmes, 1996] or [Perko, 1996]. Example 13. Let us fix d = 0.2, = 0.1, r =
0.16, β = 1.22. Then a = 0.1311, h = 0.2951, c = 0.2459, H = 0.3771 − 0.2951σ. If we vary σ from 0 to 1, we have a1 = −0.2951(σ + 1.0865)(σ + 0.3579) < 0, a0 = 0.0726σ 2 + 0.0806σ + 0.2075 > 0. Further, we have 2a0 + a1 > 0 if 0 ≤ σ < 0.7845 and 2a0 + a1 < 0 if 0.7845 < σ ≤ 1. Then it is easy to see that Ei , i = 1, 2, exist when 0.4898 < σ ≤ 1. Note that r1 = 0.0053(σ + 2.8579)(σ + 2.3494) × (σ − 0.6473)(σ − 0.9822). Thus, E2 is stable when 0.6473 < σ ≤ 1 and is unstable when 0.4898 < σ < 0.6473. This means that a Hopf bifurcation may occur when σ passes through 0.6473. Since ξ1 = −0.125 when σ = 0.6473, it follows from Theorem 12 that (16)
admits an unstable limit cycle when σ increases from 0.6473. If 0 ≤ σ < 0.4898, the model has no positive equilibrium or has one positive saddle point, and the dynamical behavior is similar to those indicated in Figs. 3 and 4. If 0.4898 < σ < 0.6473, computer simulations show that all positive orbits except two endemic equilibria and the stable manifolds of E1 tend to the infection-free equilibrium as t tends to ∞ (see Fig. 5). This means that the infection dies out although two endemic equilibria occur. When σ increases from 0.6473, we have an unstable limit cycle inside which all orbits approach E2 as t tends to ∞ and outside which all orbits tend to the infection-free equilibrium as t tends to ∞ (see Fig. 6). Further, as σ increases, the limit cycle expands and a homoclinic orbit occurs when σ = 0.667938 (see Fig. 7). Then, if we further increase σ, the homoclinic orbit is broken, and the stable manifolds of E1 split the feasible region into two parts. If an initial position lies in the lower part, the orbit tends to the infection-free equilibrium. If an initial position lies in the upper part, the orbit tends to the endemic equilibrium E2 (see Fig. 8). Hence, the infection can drive the population to extinction if the reduced reproduction coefficient σ is small, the infection dies out if the coefficient is medium, and the population coexists with parasites if the coefficient is close to 1. An unstable limit cycle
3302
W. Wang et al. 1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 5. E1 is a saddle and E2 is unstable. All the positive orbits except E1 , E2 and the stable manifolds of E1 tend to the infection-free equilibrium, where d = 0.2, = 0.1, r = 0.16, β = 1.22, σ = 0.55 and K = 2.
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
Fig. 6.
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
An unstable limit cycle occurs around E2 where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2 and σ = 0.66.
and a homoclinic orbit are the boundaries of the coexistence regions. Now, we consider two important cases to see the effects of the infection-related death and the infection-reduced birth on the dynamical behavior of (16). First, we suppose σ = 1. This means that parasites are harmless for the reproduction of the population. In this case, the boundary equilibria are the same as those discussed in last section, but positive equilibria and their Hopf bifurcations can be simplified. Indeed, we have a1 = −2H + 1, a0 = H + h − 3c, 2a0 + a1 = 2h − 6c + 1 > 0. By
Theorem 6, we see that a Hopf bifurcation can only occur when H > 1/2 and (4c − 1)H 2 + (1 − 4c)H − 2c − 6ch + h2 + 9c2 + h = 0.
(26)
Furthermore, we have ξ = a21 (2H − 1) > 0. Therefore, we can state Theorem 14. Let E2 exist. If σ = 1 and (26) holds,
then (16) has an unstable limit cycle near the bifurcation point.
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
3303
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
Fig. 7.
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
A homoclinic orbit occurs when σ = 0.667938 where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2.
1 0.9 0.8 0.7
I
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.2
0.4
0.6
0.8
1 S
1.2
1.4
1.6
1.8
2
Fig. 8. The stable manifolds of E1 split the feasible region into two regions one of which is the basin of E2 and the other one is the basin of the infection-free equilibrium, where d = 0.2, = 0.1, r = 0.16, β = 1.22, K = 2 and σ = 0.72.
Now, we consider the case where σ = 0. This case means that infected hosts do not have fecundity. Then a1 = −2c + H, a0 = 2c − H + cH, 2a0 + a1 = 2c − H + 2cH. By Theorem 6, we see that a Hopf bifurcation can only occur when H > 2c/(1 − 2c) or H < 2c where c < 1/4 is used, and 2
H − 2H + 4c = 0.
(27)
By direct calculations, we obtain ξ=− Hence, we have
H 5 (H − 2) . 8
(28)
Theorem 15. Assume that E2 exists and σ = 0.
Then (16) has an unstable limit cycle near the Hopf bifurcation point if H > 2 and has a stable limit cycle near the Hopf bifurcation point if 2c/(1 − 2c) < H < 2 or H < 2c. Example 16. Let us fix r = 1, d = 0.1, σ =
0, = 0.1, β = 0.9 in (3). Then we have a = 1.1111, h = 1.2222, c = 0.2222, H = 1.3333. It is easy to see that condition (1) in Theorem 4 holds. Since 0.5714 = 2c/(1 − 2c) < H < 2, it follows that (3) admits a stable limit cycle when β is close to 0.9 (see Fig. 9). By computer simulations, we
3304
W. Wang et al. 1.5
I
1
0.5
0 0
Fig. 9.
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
2
A stable limit cycle around E2 where r = 1, d = 0.1, = 0.1, σ = 0, K = 2 and β = 0.91.
see that the periodic orbit expands as β increases until 0.9625039, and meets a homoclinic orbit at β = 0.9625039 (see Fig. 10, which is produced by MatCont [Dhooge et al., 2004]). Then, the homoclinic orbit is broken and all positive orbits except the infected equilibria and the stable manifolds of a positive saddle tend to the disease-free equilibrium when 0.9625039 < β < 1.2. This means that the infection dies out when β lies in this interval. Hence, higher contact coefficient can eradicate the infection. If we increase β from 1.2, two characteristic directions θ1 and θ2 of (0, 0) occur and there are bistable regions. Positive orbits in the lower part of the feasible region tend to the disease-free equilibrium, while positive orbits in the upper part
1.2 1
I
0.8 0.6
of the feasible region tend to (0, 0), which means that the infection drives the population to extinction. Furthermore, it is easily verified that the disease-free equilibrium is globally stable if 0 < β < 0.8, and the infection can be persistent for some initial values if 0.8 < β < 0.9625039. Hence, the model goes through the stages of extinction of the infection, persistence of infection, the extinction of infection, and the extinction of both infection and host population, as β increases from 0 to ∞.
3.3. Higher codimension bifurcation From the discussions above, we see that E1 and E2 coincide with E ∗ when c = 1/4. In this case, it is possible to have higher codimension bifurcations. To this end, we impose tr(J) = 0 and det(J) = 0, where J is the Jacobian matrix of (16) at an endemic equilibrium, so that we have a degenerate endemic equilibrium E ∗ . Note that y ∗ = 1/2 at the degenerate equilibrium, which implies that the determinant of the matrix J is zero. Note also that the trace of the Jacobian matrix J at the equilibrium equals zero if 1 3 −2a + H + Hσ = 0. 2 2
0.4 0.2 0
1 S
(29)
Consequently, we have E ∗ = (2a/(3 + σ), 1/2) and 0
0.5
1
1.5
S Fig. 10. Limit cycles expand as β increases from 0.91 and meet a homoclinic orbit when β = 0.9625039, where r = 1, d = 0.1, = 0.1, σ = 0, K = 2.
a1 = − a0 =
(1 + σ) (−8a + 3 + σ) , 2(3 + σ)
(1 + σ)(−8a + 3 + σ) . 4(3 + σ)
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
Theorem 17. Suppose that c = 1/4 and (29) hold.
E∗
of (16) is a cusp of codiThen the equilibrium mension 2 or higher codimension. First, we translate the equilibrium E ∗ = (2a/(3+σ), 1/2) to the origin by the transformation of X = x − 2a/(3 + σ) and Y = y − 1/2. Then we obtain dX a(1 + σ) 4a2 (1 + σ) 1 σ =− X− + X2 Y − dt 3+σ (3 + σ)2 2 2 Proof.
4aσ XY + o((X, Y )2 ), 3+σ 1 σ a (1 + σ) dY = + X+ Y + σXY dt 4 4 3+σ −
+
If Z = −(a(1 + σ)X/(3 + σ)) − 4 (a2 (1 + σ)Y/ (3 + σ)2 ), it becomes 1−σ 2 σ(3 + σ) dX =Z− X + XZ dt 2 a(1 + σ)
(30)
(3 + σ) (−3 − σ + 4a + 4aσ) 2 Z 8a2 (1 + σ)
+ o((X, Z)2 ). (3 + σ) (−3 − σ + 4a + 4aσ) XZ . P =Z+ 8a2 (1 + σ) Then (30) is transformed into 1−σ 2 dX =P− X dt 2 (3 + σ) (−4a + 4aσ + 3 + σ) XP a2 (1 + σ)
+ o((X, P )2 ), 1 + σ 2 4a(1 + σ) − 3 − σ dP = X − XP dt 8 4a + o((X, P )2 ). Set u=X+
(−3 − σ + 4a + 4aσ)X 2 , 8a
dv = dt
σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a uv 8a2 (1 + σ)
(32)
1−σ 2 u + o((u, v)2 ). 2 If σ = 1 and −
then E ∗ is a cusp of codimension 2, i.e. there is a Bogdanov–Takens bifurcation [Chow et al., 1994, Sec. 4.1]. If σ = 1 or σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a = 0,
One implication of the bifurcation of codimension 3 or higher codimension is that (16) has two limit cycles and homoclinic bifurcation under suitable perturbations.
4. Discussion
Set
+
du = v + o((u, v)2 ), dt
it follows from [Chow et al., 1994, Sec. 5.3] that E ∗ is a cusp of codimension 3 or higher codimension.
+ o((X, Z)2 ),
−
(1 + σ)X 2 . 8 Then (31) becomes v=P+
σ 2 (4a + 1 + 2a2 ) + σ(6 + 4a2 + 8a) + 2a2 + 9 − 12a = 0,
4a(1 + σ) − 3 − σ 2 Y + o((X, Y )2 ). 2(3 + σ)
1 + σ 2 4a(1 + σ) − 3 − σ dZ = X − XZ dt 8 4a
3305
(31)
Host extinction and host oscillations are reported in papers [Ebert et al., 2000; Krukonis & Schaffer 1991; Kendall et al., 1999]. Ebert et al. [2000] showed that stochastic factors are possible causes for these. The paper by Hwang and Kuang [2003] indicates that a standard incidence can lead to the extinction of a host population. This means that deterministic factors are also causes for the host extinction. In this paper, we have adopted the incidence βSI 2 /N 2 to simulate the situation where the infection probability of a susceptible individual depends linearly upon the number of exposures to infected individuals. Theorem 1 shows that the host population cannot go to extinction if β < (d + r)(1 + σ) + . This means that the population is persistent if β is small. If β > (d + r)(1 + σ) + , we have shown that the model admits saddle node bifurcations, Hopf bifurcations and codimension 2 or higher codimension
3306
W. Wang et al.
bifurcations. Numerical calculations indicate that the infection can die out even though there are two positive equilibria. The more interesting result is that the model admits bistable regions of initial states for many parameter values, in which the infection dies out if initial states lie in one region and the population and parasites coexist if initial states lie in the other region, and admits a region of initial states such that the population and parasites die out together in this region for suitable parameters. Since the extinction of the population or the infection is dependent upon initial states of the population and parasites, this seems more reasonable in reality. In contrast, either the infection dies out and the population is persistent, or both parasites and the population go to extinction together, if the standard incidence is used [Hwang and Kuang, 2003]. Hence, the correlation of infection probability of a susceptible individual to multiple exposures to infectives is a cause for bistable states. Another scenario is also interesting. With the nonlinear incidence incorporating the effect of multiple exposures, stable and unstable limit cycles can occur in model (5). Since the model with the standard incidence does not have a limit cycle (see [Hwang & Kuang, 2003]), this means that the infection probability related to multiple exposures to infected individuals is one source of oscillations of the host population. If parasites are harmless for the reproduction, i.e. σ = 1, we have shown that a Hopf bifurcation is subcritical, i.e. an unstable limit cycle is produced. For σ < 1, we find that the model has stable limit cycles. It seems that the reduction coefficient σ of birth rate of infected hosts is important for producing a stable periodic oscillation. If we decrease σ from 1 to 0, we have found in Example 13 that the model undergoes the coexistence of the population and the infection in one region of initial states, extinction of the infection, and extinction of the population and the infection in one region of initial states. Hence, medium reduction of birth rate of infected hosts seems better for infection control. The higher reduction of birth rate of infected hosts is dangerous for the population, and lower reduction increases the risk of infection. We have also found in Example 16 that the model goes through the stages of extinction of the infection, persistence of the infection in one region of initial states, extinction of the infection, and extinction of both infection and the population in one
region of initial states, as β increases from 0 to ∞. These extinction and persistence switches reveal the effects of β. High β is dangerous for the population and suitably higher β can also eradicate the infection.
References Begon, M., Bennett, M., Bowers, R. G., French, N. P., Hazel, S. M. & Turner, J. [2002] “A clarification of transmission terms in host-microparasite models: Numbers, densities and areas,” Epidemiol. Infect. 129, 147–153. Berezovskaya, F. Karev, G. & Arditi, R. [2001] “Parametric analysis of the ratio-dependent predator–prey model,” J. Math. Biol. 43, 221–246. Bouma, A., De Jong, M. C. M. & Kimman, T. G. [1995] “Transmission of pseudobabies virus within pig-populations is independent of the size of the population,” Pre Vet Med. 23, 163–172. Castillo-Chavez, C. & Yakubu, A. A. [2001] “Discretetime S-I-S models with complex dynamics,” Nonlin. Anal. 47, 4753–4762. Chow, S.-N., Li, C. Z. & Wang, D., [1994] Normal Forms and Bifurcation of Planar Vector Fields (Cambridge University Press, Cambridge). De Jong, M. C. M., Diekmann, O. & Heesterbeek, J. A. P. [1995] “How does transmission depend on population size?” in Human Infectious Diseases, Epidemic Models, ed. Mollison, D. (Cambridge University Press, Cambridge), pp. 84–94. Deredec, A. & Courchamp, F. [2003] “Extinction thresholds in host-parasite dynamics,” Ann. Zool. Fennici 40, 115–130. Derrick, W. R. & van den Driessche, P. [1993] “An infection transmission model in a nonconstant population,” J. Math. Biol. 31, 495–512. Derrick, W. R. & van den Driessche, P. [2003] “Homoclinic orbits in an infection transmission model,” Discr. Contin. Dyn. Syst. Ser. B. 3, 299–309. Dhooge, A., Govaerts, W. & Kuznetsov, Y. A. [2004] Limit cycles and their bifurcations in MatCont. http://allserv.UGent.be/ajdhooge. Diekmann, O. & Heesterbeek, J. A. P. [2000] Mathematical Epidemiology of Infectious Diseases. Model Building, Analysis and Interpretation, Wiley Series in Mathematical and Computational Biology (John Wiley & Sons, Ltd., Chichester). Ebert, D., Lipsitch, M. & Mangin, K. L. [2000] “The effect of parasites on host population density and extinction: Experimental epidemiology with Daphnia and six microparasites,” Amer. Naturalist 156, 459– 477. Gao, L. & Hethcote, H. W. [1992] “Disease transmission models with density-dependent demographics,” J. Math. Biol. 30, 717–731.
Bifurcations in a Host-Parasite Model with Nonlinear Incidence
Grenfell, B. T. & Dobson, A. P. (eds.) [1995] Ecology of Disease in Natural Populations (Cambridge University Press, Cambridge). Guckenheimer, J. & Holmes, P. J. [1996] Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer-Verlag, NY). Hethcote, H. W. [2000] “The mathematics of infectious diseases,” SIAM Rev. 42, 599–653. Hwang, T. W. & Kuang, Y. [2003] “Deterministic extinction effect of parasites on host populations,” J. Math. Biol. 46, 17–30. Kendall, B. E., Briggs, C. J., Murdoch, W. W., Turchin, P., Ellner, S. P., McCauley, E., Nisbet, R. M. & Wood, S. N. [1999] “Why do populations cycle? A synthesis of statistical and mechanistic modelling approaches,” Ecology 80, 1789–1805. Krukonis, G. & Schaffer, W. M. [1991] “Population cycles in mammals and birds: Does periodicity scale with body size?” J. Theor. Biol. 148, 469–493. Liu, W. M., Levin, S. A. & Iwasa, Y. [1986] “Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models,” J. Math. Biol. 23, 187–204. Liu, W. M., Hethcote, H. W. & Levin, S. A. [1987] “Dynamical behavior of epidemiological models with nonlinear incidence rates,” J. Math. Biol. 25, 359–380. Lizana, M. & Rivero, J. [1996] “Multiparametric bifurcations for a model in epidemiology,” J. Math. Biol. 35, 21–36. May, R. M. & Anderson, R. M. [1979] “Population biology of infectious diseases. II,” Nature 280, 455–461.
3307
McCallum, H. [1995] “Modelling wildlife-parasite interactions to help plan and interpret field studies,” Wildlife Res. 22, 21–29. Nold, A. [1980] “Heterogeneity in disease-transmission modeling,” Math. Biosci. 52, 227–240. Perko, L. [1996] Differential Equations and Dynamical Systems (Springer-Verlag, NY). Polking, J. C., http://math.rice.edu/dfield/6.5. Ruan, S. & Wang, W. [2003] “Dynamical behavior of an epidemic model with a nonlinear incidence rate,” J. Diff. Eqs. 188, 135–163. van den Driessche, P. & Watmough, J. [2002] “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Math. Biosci. 180, 29–48. van Nouhuys, S. & Ehrnsten, J. [2004] “Wasp behavior that leads to uniform parasitism of a host available only a few hours per year,” Behav. Ecol., in press. Wang, W. & Ruan, S. [2004] “Simulating the SARS outbreak in Beijing with limited data,” J. Theor. Biol. 227, 369–379. Wang, W. & Zhao, X.-Q. [2004] “An epidemic model in a patchy environment,” Math. Biosci. 190, 97– 112. Xiao, D. & Ruan, S. [2001] “Global dynamics of a ratiodependent predator–prey system,” J. Math. Biol. 43, 268–290. Zhang, Z., Ding, T., Huang, W. & Dong, Z. [1991] Qualitative Theory of Differential Equations, Trans Mathematical Monographs, Vol. 101 (Amer. Math. Soc., Providence).
