Bifurcation analysis in models for vector-borne diseases with logistic growth.

We establish and study vector-borne models with logistic and exponential growth of vector and host populations, respectively. We discuss and analyses the existence and stability of equilibria. The model has backward bifurcation and may have no, one, or two positive equilibria when the basic reproduction number R 0 is less than one and one, two, or three endemic equilibria when R 0 is greater than one under different conditions. Furthermore, we prove that the disease-free equilibrium is stable if R 0 is less than 1, it is unstable otherwise. At last, by numerical simulation, we find rich dynamical behaviors in the model. By taking the natural death rate of host population as a bifurcation parameter, we find that the system may undergo a backward bifurcation, saddle-node bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation, and cusp bifurcation with the saturation parameter varying. The natural death rate of host population is a crucial parameter. If the natural death rate is higher, then the host population and the disease will die out. If it is smaller, then the host and vector population will coexist. If it is middle, the period solution will occur. Thus, with the parameter varying, the disease will spread, occur periodically, and finally become extinct.

1. Introduction

The dynamical modeling and studies of vector-borne disease have been in many literatures. The earliest studies originated from the Ross-Macdonald model for malaria [1] in 1911 and major extensions are described in Macdonald's 1957 book [2]. For malaria, there aremany literature [3-7]. Research literatures on dengue fever have [8-12], and so forth. For West Nile virus (WNV), there are also many literatures, [13-18], and so forth. Here, we classify vector population as susceptible (S ) and infected (I ) and classify the host population as susceptible (S ), infected (I ), and recovered R . In these literatures, birth function of susceptible host and vector population and incidence rate are taken as Table 1. Parameters involved in Table 1 for the dynamics of vector-borne diseases are in Table 2. For malaria, dengue fever, and West Nile virus, the three typical vector-borne diseases, the dynamics of most of the related models can be characterized by the reproduction number, but some compartmental models for these three mosquito-borne diseases undergo a backward bifurcation, an important feature of mosquito-borne diseases.

Table 1

The birth and incidence rates for different vector-borne disease models.

Disease	Literature	Birth		Incidence rate
		S h	S v	Host	Vector
	[3]	λh N h	λv N v	β_hI_vS_hN_h	β_vI_hS_vN_h
Malaria	[5, 6]	Λh + λh N h	λv N v	β_hI_vS_hσ_vN_v+σ_hN_h	β_vI_hS_vσ_vN_v+σ_hN_h
	[7]	(1 − ϕ)Λh	λv N v	β_hI_vS_hN_h	β_vI_hS_vN_h
	[8]	λh N h	Λh	β_hI_vS_hN_h+m	β_vI_hS_vN_h+m
Dengue fever	[9]	λh N h	λv N v	β_hI_vS_hN_v	β_vI_hS_vN_h
	[10]	λh N h	λv N v	β_hI_vS_hN_h+m	β_vI_hS_vN_h+m
	[13]			β_1I_vS_hN_h	β_vI_hS_vN_h
WNV	[14]	Λh	Λv	β_hI_vS_hN_h	β_vI_hS_vN_h
	[15, 18]	Λh	λv N v	β_hI_vS_hN_h	β_vI_hS_vN_h

Table 2

Parameters involved in the models for the dynamics of vector-borne diseases.

Parameters in the models	Vector	Host
Per capita birth rate	λv	λh
The fraction of the infective immigrants		ϕ
Number of times one mosquito would want to bite humans per unit time	σh
The maximum number of mosquito bites a human can have per unit time		σv
Recruitment rate per unit of time	Λv	Λh
Per capita birth rate	ψv	ψh
Vertical transmission rate in new birth of vectors	p
Natural death rate	d v	d h
Disease-induced death rate		αh
Recovery rate		γh
Biting rate (the average number of bites per mosquito per day)	b
Transmission probability (from infected vectors to host)		βh
Transmission probability (from infected host to vectors)	βv

For compartmental models for malaria, dengue fever, and WNV, different forms of birth functions have been used; we summarize and find that the logistic birth function which describes the self-limiting growth of a host population was rarely considered. In the paper, we will build the following model and analyse the dynamic behaviors: In the model (1), the host and vector populations satisfy the following equations, respectively: It means that the total number of the vector populations is a constant and denotes N = N . Furthermore, using the relations S + I + R = N and S = N − I , we lead to the following four-dimensional nonlinear system of ODEs:

One can verify that the positive cone is positive invariant.

The remainder of this paper is organized as follows. In Section 2, we consider the number and existence of equilibria. Section 3 studies local stability of the disease-free equilibrium and gives some examples to explain the existence of endemic equilibria and bifurcation with the basic reproduction number changing. The paper concludes with a brief discussion in Section 4.

2. Existence of Equilibria

To study the existence of equilibria, let the right side of each of the four differential equations be equal to zero in system (3), obtaining the following equations:

From the above equations, one can verify that system (3) has a disease-free equilibrium at E 0 : (S , I , N , I ) = (K(1 − d /r), 0, K(1 − d /r), 0). The linear stability of E 0 is governed by the basic reproduction number R 0 which can be found from the next generation matrix [19] for the system. Noting that the model has infected populations, namely, I and I , it follows that, using the notation of van den Driessche and Watmough [20], the matrices F and V, for the new infection terms and the remaining transfer terms, respectively, are given in partitioned form by where Here, F is a nonnegative matrix of rank 2 and V is nonsingular. It is easy to show that, at steady state, the spectral radius (dominant eigenvalue) of the nonnegative matrix FV −1, denoted by ρ(FV −1), is equal to ρ(F 1 V 1 −1); hence, One can see that R 0 is the geometric mean of the number of new diseased generated in both the host and vector population, respectively. The first part is the number of host infections caused by one infected vector; the second part is the number of vector infections caused by one infected host. Discussion in detail on basic reproductive number of host-vector disease models see literature [16].

In order to find positive equilibria, it follows from the first and last two equations of (5), we have Substituting (10) into the second equation of (5), that is, if an endemic equilibrium exists, its N coordinate satisfies where

N should satisfy 0 < N < N = K(1 − d /r) to ensure a positive equilibrium. Then we will consider the existence of positive roots for (11) from the following two cases.

2.1. Case 1: a 1 ≠ 0, a 2 = 0

In this case, system (3) for vector-borne diseases is reduced to For system (13), N satisfies where

Now we consider the case 0 ≤ p < 1. One can see that it is the number of roots in the interval (0, K(1 − d /r)) that determines the number of positive equilibria for model (13). Note If R 0 < 1, F 1(K(1 − d /r)) < 0, F 1(N ) = 0 has a unique root in (0, K(1 − d /r)). If R 0 > 1, F 1(K(1 − d /r)) > 0, then there may be two positive equilibria or no positive equilibrium. We have the following theorem.

Theorem 1

If R 0 < 1, system (13) has a unique positive equilibrium E*(S *, I *, N *, I *).

Theorem 2

If R 0 > 1, one has the following.

If −rB 1 ≥ K(r − d )B 2, then system (13) has no endemic equilibrium.

If −rB 1 < K(r − d )B 2, then

when B 1 2 < 4B 2 B 0, system (13) has no endemic equilibrium;

when B 1 2 = 4B 2 B 0, system (13) has a unique repeated endemic equilibrium E** = (S **, I **, N **, I **);

when B 1 2 > 4B 2 B 0, system (13) has two endemic equilibria E 1*(S ∗1, I ∗1, N ∗1, I ∗1), E 2*(S ∗2, I ∗2, N ∗2, I ∗2), where N ∗1 < N ** < N ∗2.

2.2. Case 2: a 1 ≠ 0, a 2 ≠ 0

For F(N ) of system (3), If R 0 < 1, F(K(1 − d /r)) < 0 and R 0 > 1, F(K(1 − d /r)) > 0. Then we will analyze the existence of positive roots of F(N ) = 0. The derivative of F(N ) is Denote Then we have the following theorem.

Theorem 3

If R 0 < 1, one has the following.

System (3) has no endemic equilibrium, if any one of the following conditions holds:

A 2 > 0 and A 1 < 0,

A 2 < 0, A 1 > 0, and K(1 − d /r) < N ,

A 2 < 0, A 1 > 0, K(1 − d /r) > N , and F(N ) > 0.

System (3) has two endemic equilibria, if A 2 < 0, A 1 > 0, K(1 − d /r) > N , and F(N ) < 0.

Theorem 4

If R 0 > 1, one has the following.

System (3) has a unique endemic equilibrium, if any one of the following conditions holds:

A 2 > 0 and A 1 < 0;

A 2 < 0, A 1 > 0, A 2 2 ≤ 3A 3 A 1;

A 2 < 0, A 1 > 0, A 2 2 > 3A 3 A 1, and K(1 − d /r) < N ;

A 2 < 0, A 1 > 0, A 2 2 > 3A 3 A 1, K(1 − d /r) > N , F(N ) < 0;

A 2 < 0, A 1 > 0, A 2 2 > 3A 3 A 1, K(1 − d /r) > N , F(N ) > 0, and F(N ) > 0.

System (3) has three endemic equilibria, if A 2 < 0, A 1 > 0, A 2 2 > 3A 3 A 1, K(1 − d /r) > N , F(N ) > 0, and F(N ) < 0.

It is easy to prove Theorem 4, we will omit the proof. In the following part, we will give two examples to show changes of equilibria with the basic reproduction number R 0 by numerical simulation when a 2 = 0 and a 2 ≠ 0.

Example 5

Fix parameters r = 0.8, K = 10, a 1 = 1, a 2 = 0, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2. Figure 1 shows that there is always a unique endemic equilibrium when R 0 < 1. As R 0 increases, the equilibrium will be keep. If R 0 equals one, a new positive equilibrium is branched from the disease-free equilibrium. If R 0 continue to be increased, the two positive equilibria may be repeated and then disappeared if R 0 exceeds a certain value, that is, there are two positive equilibria or no positive equilibrium when R 0 > 1 if a 2 = 0.

Figure 1

Bifurcation curves in (R 0, I ) plane when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Example 6

Take r, K, a 1, d , γ , α , p, d , β 2, and N as same values of Example 5 and a 2 = 0.1. From Figure 2, we find that the number of positive equilibria is from zero to two with R 0 increasing when R 0 < 1. If R 0 is more than one, it is found that the two positive equilibria still be keep and there is a new endemic one to be branched from the disease-free equilibrium. Now system (3) has three endemic equilibria. Continuing to increase R 0, we find that the larger two equilibria will be repeated and then disappear and the smaller equilibrium always be keep. Let a 2 = 0.3; from Figure 3, it can be obtained that there is no positive equilibrium when R 0 < 1. With R 0 increasing, a positive equilibrium is branched from the disease-free one; then three ones appear and the larger two ones disappear and the smaller one lasts.

Figure 2

Bifurcation curves in (R 0, I ) plane when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.1, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Figure 3

Bifurcation curves in (R 0, I ) plane when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

3. Stability and Bifurcation Analysis

Then we study the stability of the steady states of system (3) and possible bifurcations. Firstly, we discuss the stability of disease-free equilibrium E 0. At this equilibrium point, the Jacobian matrix is where The characteristic polynomial of (20) is It means that or

The roots of the two quadratic polynomial have negative real parts if and only if their coefficients are positive; that is, R 0 < 1. Therefore, the disease-free equilibrium E 0 is always locally asymptotically stable when R 0 < 1, and it is unstable if R 0 > 1.

For any positive endemic equilibrium E , the Jacobian matrix of (3) becomes where Here, it is difficult to determine a sign for the real part of eigenvalue of J (E ), so we will give the dynamical behaviors for system (3) by numerical simulation.

Example 7

Fix parameters r = 0.8, K = 10, a 1 = 1, a 2 = 0, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2. In Figure 1, we find that a unique positive equilibrium is unstable when R 0 < 1. With R 0 increasing and exceeding one, the unstable positive equilibrium still be keep and a new positive one branching from the disease-free one is stable. After R 0 is more than a threshold, the two equilibria will disappear, and all of trajectories tend to equilibrium E(S , I , R , S , I ) = (0,0, 0, N , 0), which means that the host population will be extinct If we fix R 0 and let d change, we find that Hopf bifurcation maybe occurs from the blue curve of Figure 4 for system (13).

Figure 4

Bifurcation curves in (d , I ) plane when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0,0.1,0.3, d = 0.036, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Example 8

Fix parameters r = 0.8, K = 10, a 1 = 1, β = 1, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2. Then we draw the bifurcation figure with d changed when a 2 = 0,0.1,0.3, respectively. We use blue, black, and brown curves to represent the bifurcation curves when a 2 = 0,0.1,0.3. In Figure 4, we firstly choose d as bifurcation parameter and draw the bifurcation curve (blue (a 2 = 0), black (a 2 = 0.1), and brown curves (a 2 = 0.3)) I with d . When a 2 = 0, we find that R 0 is always more than 1 no matter what the value of d is, which means that system (3) has no or two positive equilibria. From the blue curve of Figure 4, we can see that there are a limit point (LP) at d = 0.208198 and a Hopf point (H) at d = 0.012541 whose first Lyapunov coefficient is 4.560482e − 001. When a 2 = 0.1, it is found that two limit points at d = 0.235734 and 0.005803 and a Hopf point at d = 0.016349 whose first Lyapunov coefficient is 6.142143e − 001 for system (3). This is consistent with the existence of the positive equilibrium. If we take a 2 = 0.3, by numerical simulation, we see that there are two limit points at d = 0.296694 and 0.041986, two neutral saddle points at d = 0.236797 and 0.177735, and two Hopf points at d = 0.043395 (first Lyapunov coefficient 5.036862e + 000) and d = 0.748616 (first Lyapunov coefficient −1.827258e − 001). Then take a limit point as initial point and d and a 2 as bifurcation parameters and obtain a fold curve (red curve in Figure 4). In the fold curve, there are a cusp bifurcation point at d = 0.434317, a 2 = 0.631081, and two Bogdanov-Takens bifurcation points at d = 0.414918, a 2 = 0.608913 and d = 0.074370, a 2 = 0.356609. At the cusp bifurcation point two branches of saddle-node bifurcation curve meet tangentially, forming a semicubic parabola. For nearby parameter values, the system has three equilibria which collide and disappear pairwise via the saddle-node bifurcations. For the Bogdanov-Takens (BT) bifurcation nearby parameter values, the system has two equilibria which collide and disappear via a saddle-node bifurcation. In the end, take any one of Bogdanov-Takens bifurcation points as initial point and draw Hopf bifurcation curve (the green curve of Figure 4). It shows that there is a Zero-Neutral Saddle (ZH) at d = 0.085991 and a 2 = 0.370274. Thus, we can obtain that all of points being saddle between two Bogdanov-Takens bifurcation points and Hopf points locate on both sides of Bogdanov-Takens bifurcation points from Hopf bifurcation curve of Figure 4.

In the following example, we will analyse how to change the state variables I and I with time t.

Example 9

Take r, K, β 1, a 1 γ ,α , p,d ,β 2, and N as the same with the above example and a 2 = 0.3. We see how to change for trajectories I and I with time t when d = 0.2,0.3, 0.6,0.7486, and 0.8, respectively. From Figures 5 and 9, we can find that all of trajectories tend to an endemic equilibrium when d is smaller, and all of trajectories tend to an equilibrium (S , I , R , S , I ) = (0,0, 0, N , 0) when d is larger, that is, the host population will be distinct, and the vector population will be keep. If d = 0.3,0.6,0.7486, respectively, from Figures 6, 7, and 8, we find that system (3) has period solution. Furthermore, from these figures, we can see that the numbers of the infective host population and the infective vector population are on the decline with the natural death rate increase, but the rate of the decline for host population is much faster than that for vector population. Finally, the host population decays to zero.

Figure 5

Trajectory picture with time t when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.2, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Figure 9

Trajectory picture with time t when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.8, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Figure 6

Trajectory picture with time t when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.3, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Figure 7

Trajectory picture with time t when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.6, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

Figure 8

Trajectory picture with time t when r = 0.8, K = 10, β 1 = 1, a 1 = 1, a 2 = 0.3, d = 0.7486, γ = 0.06, α = 0.36, p = 0.007, d = 0.06, β 2 = 0.26, and N = 2.

4. Discussion

In previous literatures, the dynamic behaviors of models on vector-borne disease may have backward bifurcation and are not complex. In this paper, we consider the behaviors of a model with standard incidence rate and the logistic birth function for host population and the exponential birth rate for vector population. The rich dynamic behaviors are found. We find that the complex behaviors will not occur if the birth function is not logistic. Taking the natural death rate of host population as a bifurcation parameter, we find that the system may undergo a backward bifurcation, saddle-node bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation, and cusp bifurcation with the saturation parameter varying. We find that the natural death rate still is a key factor. If the natural death rate is higher, then the host population and the disease will die out; if it is smaller, then the host and vector population will coexist; that is, the disease will spread. If it is middle, the period solution will occur. In addition, the existence of equilibria of our model is different from previous model. If the incidence rate function is not affected by the number of vector populations, we find that there is always a unique positive equilibrium when the basic reproduction number R 0 is less than one, and there are two endemic equilibria or no positive equilibrium when R 0 is more than one. If the incidence rate function depends on the number of vector populations, we obtained that there may be two positive equilibria or no positive equilibrium when R 0 < 1 and that there exist one or three endemic equilibria when R 0 > 1 where the unique endemic one is always unstable.