The impact of predators of mosquito larvae on Wolbachia spreading dynamics

Dengue fever creates more than 390 million cases worldwide yearly. The most effective way to deal with this mosquito-borne disease is to control the vectors. In this work we consider two weapons, the endosymbiotic bacteria Wolbachia and predators of mosquito larvae, for combating the disease. As Wolbachia-infected mosquitoes are less able to transmit dengue virus, releasing infected mosquitoes to invade wild mosquito populations helps to reduce dengue transmission. Besides this measure, the introduction of predators of mosquito larvae can control mosquito population. To evaluate the impact of the predators on Wolbachia spreading dynamics, we develop a stage-structured five-dimensional model, which links the predator-prey dynamics with the Wolbachia spreading. By comparatively analysing the dynamics of the models without and with predators, we observe that the introduction of the predators augments the number of coexistence equilibria and impedes Wolbachia spreading. Some numerical simulations are presented to support and expand our theoretical results.


Introduction
Dengue fever, the fastest spreading mosquito-borne disease in tropics and subtropics, causes 390 million confirmed cases annually in the world, with about 36 thousand deaths [7].The sterile insect technique (SIT) [40] and the incompatible insect technique (IIT) [10] are two effective and eco-friendly weapons for stemming the spread of dengue fever, which respectively employ the radiation-treated or the endosymbiotic bacterium Wolbachia-infected male mosquitoes (we refer to them as sterile mosquitoes hereafter) reared in the labs or mosquito factories to sterilize wild mosquitoes, and hence to suppress the local mosquito population.After the sterile mosquitoes being released into the habitat of the wild mosquitoes, the wild and sterile mosquitoes will show interactive behaviours.To explore the interactive dynamics between wild and sterile mosquitoes and hence to help the staff working in the line of SIT or IIT to design cost-effective strategies for releasing sterile mosquitoes, a wealth of mathematical models with and without time delay have been developed and investigated, we refer to [4, 15-17, 21, 36-39, 43-45] and references cited therein.
Mosquitoes undergo complete metamorphosis going through four distinct phases of development during their whole-life: eggs, larvae, pupae and adults, and the first three stages (called larvae) are in water, while the last stage is in the air.Accordingly, when applying SIT or IIT for suppressing wild mosquitoes and preventing mosquito-borne diseases, it is rational to take into account the effects of the stage structures of wild mosquitoes.To assess the impact of the stage structures of wild mosquitoes on their suppression dynamics, Ai et al. [1] proposed the following model Here, J = J(t), A = A(t) denote the numbers of mosquito larvae and adult mosquitoes at time t, respectively, and all the parameters are positive and their implications are referred to Table 1.Through analysing the model, the authors perceived that the stage-structured two-dimensional system can exhibit much more complicated dynamics than those onedimensional models.Nevertheless, both SIT and IIT require successive releases of sterile mosquitoes for obtaining the sustained suppression effect, which consume a lot of manpower, mosquitoes and financial resources.On the other hand, experiment data show that some strain of Wolbachia can inhibit the replication of the virus in mosquitoes [2,33], and induce the mechanisms of maternal transmission and cytoplasmic incompatibility (CI) [31,34], which leads to an economical alternative strategy named mosquito population replacement: to release moderate Wolbachia-infected mosquitoes (pure females or mixtures of females and males) to invade wild mosquito population and reduce the spread of mosquito-borne diseases.Once Wolbachia-infected mosquitoes have been released into the target area, Wolbachia start to invade the local mosquito population.Deeper insights into the spreading dynamics contribute a lot to the stable establishments of Wolbachia.To that effect, a battery of mathematical models have been formulated and analysed, see, for example, [14,25,27,35,41,42] and references cited therein.
To develop a model characterizing the spreading dynamics of Wolbachia in mosquitoes while mathematically tractable, we suppose that the following assumptions hold hereafter.Based on the above model and assumptions, we arrive at the following system: which depicts the Wolbachia spreading dynamics in mosquito population.Here, J U = J U (t), J I = J I (t), A U = A U (t) and A I = A I (t) are the numbers of uninfected larvae, infected larvae, uninfected adults and infected adults at time t, respectively.Furthermore, introducing predators (such as the arroyo chub) of mosquito larvae is also a feasible solution to control mosquito population density.This method is compatible with the aquatic environments, and there have been a body of successful cases [6,11,32].Since mosquito larvae can't fly and the space of their habitats are limited, introducing aquatic predators to control the density of mosquito larvae is reasonable and economical [26,46].
In this work we focus on the combination of the above two strategies for preventing the prevalence of dengue fever, that is, release a specific quantity of predators of mosquito larvae to a place where Wolbachia-infected and uninfected mosquitoes coexist.Let P = P(t) represent the number of predators of mosquito larvae at time t, and assume that the predators will prey on infected and uninfected mosquito larvae equally without any obvious distinguish.Considering the Wolbachia spreading dynamics in mosquitoes and the interactive dynamics of predators and mosquito larvae with Michaelis-Menten (or Holling-II) functional response [8,12,13,19,23], we develop the following system: where the explanations of the positive parameters c, a, m and δ are also given in Table 1.
In this study, we are dedicated to figuring out the impact of the predators of mosquito larvae on Wolbachia spreading dynamics.In particular, we want to know whether the infection frequency of Wolbachia in mosquitoes can be improved by the introduced predators of mosquito larvae.The paper is arranged as follows.Section 2 deals with the uniqueness, positivity and boundedness of the solution of system (2).Section 3 gives the main results.More precisely, by employing the centre manifold theory, the Descartes's rule of signs and the Routh-Hurwitz criterion, the existence and stabilities of equilibria of the systems without and with predators are investigated in order.Several numerical examples are also offered to support and illustrate our theoretical results.Finally, some summaries and discussions on the systems without and with predators are provided in Section 4.
Proof: By [28,29], for any initial data (J U (0), J I (0), A U (0), A I (0), P(0)) ∈ D, system (2) has a unique and continuous solution (J U (t), J I (t), A U (t), A I (t), P(t)) defined in [0, +∞).Now, we show the positivity of the solution by contradiction.First, the fifth equation of system (2) tells us that P(t) > 0 for all t > 0. Subsequently, we consider the following four cases to prove the positivities of J U (t), J I (t), A U (t) and A I (t).
(1) Assume that both J U (t) and J I (t) have no zeros in (0, +∞), that is, J U (t) > 0 and J I (t) > 0 for all t > 0. Then from the third and fourth equations of system (2), we know that A U (t) > 0 and A I (t) > 0 for all t > 0. (2) Assume that J U (t) has at least one zero in (0, +∞) and J I (t) has no zeros in (0, +∞).
Then there exists t U > 0 such that J U (t) > 0, t ∈ (0, t U ), and J U (t U ) = 0, which implies J U (t U ) ≤ 0. Thus, the third and fourth equations of system (2) tell us that A U (t) > 0 for t ∈ (0, t U ) and A I (t) > 0 for all t > 0. In the meantime, from the first equation of system (2), we obtain (s) ds ≥ J U (0) > 0, a contradiction to J U (t U ) = 0.
From the third and fourth equations of system (2), we know that A U (t) > 0, t ∈ (0, t U ), and A U (t U ) = 0, and A I (t) > 0, t ∈ (0, t I ), and A I (t I ) = 0.
By applying the method used in case (2), we can also get the same contradiction.(4) Assume that J U (t) has no zeros in (0, +∞) and J I (t) has at least one zero in (0, +∞).
Then there exists t I > 0 such that J I (t) > 0, t ∈ (0, t I ), and J I (t I ) = 0, which shows J I (t I ) ≤ 0. Thus, the third and fourth equations of system (2) tell us that A U (t) > 0 for all t > 0 and A I (t) > 0 for t ∈ (0, t I ).In the meantime, from the second equation of system (2), we obtain This completes the proof of the positivity of the solution of system (2).Then we show that the solution (J U (t), J I (t), A U (t), A I (t), P(t)) is bounded.We first illustrate the boundednesses of J U (t) and A U (t).Since the solution is positive, from the first equation of (2), we have which implies, along with the third equation of (2), where . Thus, if JU ≤ 0, then (4) tells us that ≤ 0, which establishes the boundednesses of J U (t) and A U (t).Otherwise, (4) signals lim sup t→+∞ J U (t) ≤ JU , this, combined with the third equation of ( 2), we can also derive the boundednesses of J U (t) and A U (t).
Similarly, the boundednesses of J I (t) and A I (t) can be obtained by considering the second and fourth equations of (2).
Moreover, from the first, second and fifth equations of system (2), we know that Then the boundednesses of A U (t) and A I (t) indicate the boundedness of P(t).Thus, the proof of the boundedness of the solution of system (2) is finished.

Main results
In this section, we study the dynamics of systems ( 1) and ( 2) in order, and system (2) will be paid more attentions on, since it includes the interventions of both infected mosquitoes and predators of mosquito larvae for controlling dengue fever, and may result in much more complicated dynamics.Recall that our goal is to figure out the impact of the introduced predators of mosquito larvae on Wolbachia spreading dynamics.As a starting point, we consider the dynamics of system (1) first.

Dynamics of the system without predators
In this subsection, we investigate the dynamics of system (1), including the number of equilibria and their corresponding stabilities.

Structure of equilibria of system (1)
Obviously, owing to the remediation of (3), the origin Ē0 = (0, 0, 0, 0) is always an equilibrium of (1).For the existence of the positive equilibria, we have the following lemma.
Proof: Assume, for the sake of contradiction, that system (1) has a positive equilibrium, denoted by or, equivalently, J * U ( J * U + J * I ) = 1, a contradiction to the positivities of J * U and J * I , which establishes the proof.

Remark 3.1:
We point out here that the biological interpretation for the non-existence of positive equilibria of system (1) is that, at the potential positive equilibria, the per capita growth rate of uninfected mosquito larvae is strictly less than that of infected mosquito larvae, both of which should be equal to zero.Subsequently, we discuss the existence of the boundary equilibria of (1).First, from the second and fourth equations of (1), we know that (1) admits at most two boundary equilibria Ē1 = (0, JI , 0, ĀI ) and Ē2 = ( JU , 0, ĀU , 0), where To discuss the exact number of boundary equilibria of (1), we define and consider two cases 0 < β < β * and β > β * , for, in biology, compared with the previous two cases, the occurrence of the critical case β = β * is rare.It follows from ( 7) that Ē0 is the unique equilibrium of system (1) when 0 < β < β * , and system (1) has three equilibria: Ē0 , Ē1 and Ē2 when β > β * .Now, we turn to analyse the stabilities of equilibria of system (1).

Stabilities of equilibria of system (1)
Similar to the analyses of the existence of equilibria of system (1), we investigate the stabilities of equilibria for two cases 0 < β < β * and β > β * separately.First, if 0 < β < β * , then by (4), we have Similarly, from the proof of the boundednesses of J I (t) and A I (t), we can obtain lim t→+∞ J I (t) = 0 and lim Thus, the unique equilibrium Ē0 is globally attractive, which manifests that the wild mosquito population will die out eventually regardless of its initial population size, provided the birth rate β is less than β * .Here we provide a numerical example below to demonstrate the above results as in [24,33].which make no sense in biology and manifest the non-existence of boundary equilibria of system (1).Thus, Ē0 is the unique equilibrium of system (1) and Figure 1 shows its global attractivity.
In the following, we consider the case β > β * .Then system (1) has three equilibria: Ē0 , Ē1 and Ē2 .We explore their stabilities one by one.
(i) It is easy to see that, on the direction A I = α μ J I , the limit of the per capita growth rate of infected mosquito larvae satisfies lim which signals the instability of Ē0 .(ii) The Jacobian matrix of (1) is given by Let the parameters α, d 0 , μ and d 1 be specified as in (10).Then we get β * ≈ 0.1665.By selecting β = 0.1 ∈ (0, β * ), and plotting the graphs of J U (t), J I (t), A U (t) and A I (t) in panels (A), (B), (C) and (D), respectively, we obtain the above figure, which illustrates the global attractivity of Ē0 .

Lemma 3.2 ([5]): Consider a general system of ordinary differential equations with a parameter
Without loss of generality, it is assumed that 0 is an equilibrium of system (14) for all values of the parameter φ, that is f (0, φ) ≡ 0 for all φ.
Moreover, assume ∂x j (0, 0) is the linearization matrix of system (14) around the equilibrium 0 with φ evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts; (H2) Matrix A has a non-negative right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.
Let f k be the kth component of f and Then the local dynamics of (14) around 0 are totally determined by a and b as follows: (i) a > 0, b > 0. When φ < 0 with |φ| 1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when 0 < φ 1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium; (ii) a < 0, b < 0. When φ < 0 with |φ| 1, 0 is unstable; when 0 < φ 1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; (iii) a > 0, b < 0. When φ < 0 with |φ| 1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0 < φ 1, 0 is stable, and a positive unstable equilibrium appears; (iv) a < 0, b > 0. When φ changes from negative to positive, 0 changes its stability from stable to unstable.Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.
To apply Lemma 3.2 for the determination of the stability of Ē2 , we let X = [x 1 , x 2 , x 3 , x 4 ] T with x 1 = J U , x 2 = J I , x 3 = A U and x 4 = A I .Then (1) can be transformed as The Jacobian derivative of system (15) with β = β * is Obviously, E 2 = (x 1 , 0, x3 , 0) = ( JU , 0, ĀU , 0) is an equilibrium of system (15), and J β * (E 2 ) has a simple zero eigenvalue and all the other eigenvalues have negative real parts.Therefore, we can employ the centre manifold theory to investigate the dynamics of the system (15) near β = β * .By utilizing the method applied in [5], we can find that J β * (E 2 ) has a right eigenvector (corresponding to the zero eigenvalue), denoted by w = [w 1 , w 2 , w 3 , w 4 ] T , where and a left eigenvector (corresponding to the zero eigenvalue), given by v where It then follows from Lemma 3.2 that the dynamics of the system (15) near β = β * is fully governed by the signs of a and b.Subsequently, we estimate the values of a and b.Note that and all the other second-order derivatives are zero.The above results, together with the facts that v 1 = v 3 = 0, offer Moreover, we have Since β − β * > 0 and β is sufficiently close to β * , we obtain 0 < φ 1.Then the statement i of Lemma 15 shows that Ē2 is unstable.
We also give a numerical example below to illustrate the above results for β > β * .
To conclude, we have the following theorem.

Dynamics of the system with predators
In the previous subsection, we investigated the dynamics of system (1).In this subsection, we explore the dynamics of system (2), including the number of equilibria and their corresponding stabilities.Similar to Section 3.1, we first focus on counting the number of equilibria, and then discuss their corresponding stabilities.

Structure of equilibria of system (2)
To make the discussion on the existence of equilibria more tractable mathematically, we suppose the following two inequalities hold.
Then E * is a positive equilibrium of system (2) if and only if It is easy to see that R 6 > 0 and R 0 > 0, so we have F(x) → +∞ as x → +∞ and F(0) > 0, respectively.Moreover, direct calculations give and Then the continuity of F(x) implies that F(x) = 0 has at least two positive real roots in the interval (0, +∞), and one of which is in the interval aδ (2m − δ) , aδ (m − δ) , and the other of which is in the interval aδ (m − δ) , +∞ .This manifests that system (2) has at least one positive equilibrium.
To obtain the exact number of positive equilibria of system (2), we also need to know the signs of R 5 , R 4 , R 3 , R 2 and R 1 , according to the Descartes's rule of signs [18].Deeper computations yield That is, we have Thus, the number of changes of the signs of adjacent non-zero coefficients of F(x) is exactly two.Then the Descartes's rule of signs says that F(x) = 0 either has zero or two positive real roots, which, combined with the fact that F(x) = 0 has at least two positive real roots, imply that F(x) = 0 has exactly two positive real roots, and one of which is in the interval aδ (2m − δ) , aδ (m − δ) , and the other of which is in the interval aδ (m − δ) , +∞ .See Figure 3 for illustration.Then system (2) admits exactly one positive equilibrium with its first component lies in the interval aδ (2m − δ) , aδ (m − δ) .

Remark 3.2:
It should be stressed here that the assumptions (A1) and (A2) are only sufficient conditions to ensure system (2) to admit a unique positive equilibrium, for, the same result may be obtained under the cases when the assumption (A1) or (A2) is violated.However, for those cases, the exact numbers of positive equilibria of system (2) are too complicated to be investigated in the present time, we intend to leave them for future works.We provide a numerical example below to support the above result for the existence of a unique positive equilibrium of system (2) under the assumptions (A1) and (A2).

Stabilities of equilibria of system (2)
To investigate the stabilities of equilibria, apart from the assumptions (A1) and (A2), we suppose that the following inequality also holds.
To begin with, we consider the stability of E 0 .From the second equation of (2), we can see that, on the direction A I = α μ J I , the limit of the per capita growth rate of infected mosquito larvae satisfies lim (J U ,J I ,A U ,A I ,P)→(0,0,0,0,0) which signals the instability of E 0 .

Remark 3.3:
The instability of E 0 , together with the non-negativity of the solution, implies that there exists at least one positive component of the solution, which biologically means that the three populations can not extinct simultaneously.
Next, we check the stabilities of E i and E * , where i = 1, 2, 3, 4. The Jacobian derivative of system (2) is the matrix It follows from the assumptions (A1) and (A2) that E 1 and E 2 are unstable.
In the following, we focus on determining the stabilities of E 3 , E 4 and E * .Note that where Simple calculations, together with (16), tell us that λ 31 = −μ, λ 32 = A 1 (< 0) are two eigenvalues of J( E 3 ), and the other eigenvalues satisfy Since together with the facts that the assumptions (A1)-(A3) hold, we have It follows from the Routh-Hurwitz criterion [3] that E 3 is asymptotically stable.For the stability of E 4 , we first have , where By using ( 16) again, we find that the eigenvalues of J( E 4 ) satisfy the following equation Then (16) and the expression for B 2 imply that the two roots of are all negative.Thus, to determine the stability of E 4 , we need to know whether all roots of λ 3 + S 1 λ 2 + S 2 λ + S 3 = 0 have negative real parts.As the assumptions (A1)-(A3) hold, and similarly, by utilizing the Routh-Hurwitz criterion, we find that the equilibrium E 4 is asymptotically stable.
In the following, we discuss the stability of E * .Note that , where By letting |λI − J( E * )| = 0, we can obtain the following characteristic equation where and with It then follows from the Routh-Hurwitz criterion that the equilibrium E * is asymptotically stable if and only if the above coefficients M 1 , M 2 , M 3 , M 4 , M 5 and the relevant Routh-Hurwitz determinants (for a quintic)

Summary and discussion
Dengue fever results in considerable public health concern worldwide these years as the effects of traditional mosquito control methods are weakening [30].The endosymbiotic bacteria Wolbachia brings us new techniques for mosquito population control from two aspects: (i) the eggs from mating between uninfected females and infected males cannot hatch; (ii) Wolbachia-infected mosquitoes are less able to transmit dengue virus.The first mechanism inspired us to release only infected male mosquitoes for population suppression, and the second one prompted us to release both infected males and females for rapid population replacement.Since population suppression requires continuous release of infected male mosquitoes and consumes a lot of manpower and mosquito resources, we focused on the second measure in this study.Theoretically, one release based on rigorous mathematical calculation can lead to successful replacement.However, single control method is easy to be affected by uncertain factors, and also requires a relatively long waiting time.Therefore, we considered comprehensive strategies, and proposed a stage-structured five-dimensional model (2) in this work, which includes the Wolbachia spreading dynamics in mosquitoes and the predator-prey interaction dynamics between mosquito larvae and their predators such as the arroyo chub [6,9,11,22], aimed at evaluating the impact of the introduction of the predators of mosquito larvae on Wolbachia spreading dynamics.We started our analyses with the special system (1) which does not include predators.This model is easier to be dealt with mathematically than system (2), and the discussion on system (1) can provide research ideas for the study of system (2).We first proved that system (1) has no positive equilibria.Then under two cases 0 < β < β * and β > β * , we derived the numbers of equilibria and their corresponding stabilities.More precisely, for the case with 0 < β < β * , we found that system (1) admits Ē0 as its unique equilibrium, and it is globally attractive.While for the case with β > β * , we observed that the system has exactly three equilibria, that is, Ē0 , Ē1 and Ē2 .Moreover, the equilibria Ē0 and Ē2 are unstable, while the equilibrium Ē1 is asymptotically stable.
Subsequently, we explored the dynamics of the five dimensional system (2) under the assumptions (A1)-(A3), and observed that (2) admits five boundary equilibria, that is, the origin E 0 , and four boundary equilibria E 1 , E 2 , E 3 , E 4 .Moreover, to explore the number of positive equilibria of system (2), we first noticed that a potential positive equilibrium of (2) must satisfy (17), and then we equivalently transformed the problem of the number of positive equilibria of (2) to the problem of the number of positive real roots of (18) in a given range.By applying the Descartes's rule of signs, we finally proved that system (2) has a unique positive equilibrium E * .For the stabilities, we first found that the equilibria E 0 , E 1 and E 2 are unstable, and then by employing the Routh-Hurwitz criterion, we showed that the equilibria E 3 and E 4 are asymptotically stable, while the stability of E * is not clear in the present time.We mention here that Figure 4 just depicts the asymptotic stability of E 3 .For the asymptotic stability of E 4 , a similar figure as Figure 4 can be obtained.
Biologically, if we consider the reduced three-dimensional system composed of J U , J I and P, the equilibria E 3 , E 4 and E * describe three types of coexistence phenomena: the coexistence of infected mosquito larvae and their predators, the coexistence of uninfected mosquito larvae and their predators and the coexistence of the three of them, respectively.Compared the number of equilibria of (2) with that of (1), it is expected that the introduction of the predators will augment the number of coexistence equilibria.For system (1), Ē1 is asymptotically stable and Ē2 is unstable, which indicates that population replacement is easy to achieve.However, for system (2), both E 3 and E 4 are asymptotically stable, implying that Wolbachia-infected mosquitoes no longer have huge advantage in competition with uninfected mosquitoes.We conclude that the introduction of the predators of mosquito larvae may hinder the spread of Wolbachia in mosquitoes.The reason for the difference may be the introduced predators prey on uninfected and infected mosquito larvae equally, which declines the infection frequency of wild mosquito population and not conducive to the spread of Wolbachia in mosquitoes.Although the introduction of predators transforms E 4 from unstable equilibrium to locally stable equilibrium, it quickly reduce the number of uninfected adult mosquitoes (see Figure 5) and then reduce the risk of mosquito-borne diseases transmission.The combination of the two strategies can more effectively control mosquito populations, we will carefully study this effective combination in our future work.

( a )
Perfect maternal transmission: the offspring of infected female mosquitoes are all infected.(b) Complete CI: all the zygotic from the crossing of infected male mosquitoes with uninfected female mosquitoes will die.(c) No fitness costs: Wolbachia alters neither the birth rates of infected mosquito larvae nor the death rates of infected adult mosquitoes.(d) Gender parity: both the uninfected and infected adult mosquitoes are equally distributed in sex.

Figure 3 .
Figure 3. Let the relevant parameters in (2) are specified as in Example 3.3.Then the graph of F against x is shown in the above figure, where the left and right red dashed straight lines in the thumbnail represent x = aδ 2m−δ and x = aδ m−δ , respectively, and the intersection point of F and x-axis, marked with asterisk, corresponds to the first component of the unique equilibrium of (2) in Example 3.3.

Figure 5 .
Figure 5.We choose the parameters in Figures 2C and 4C and set A U (0) = 25.Then obtain the numbers of uninfected adult mosquitoes with and without predators.

Table 1 .
Model parameters and their interpretations.
M 3 , M 4 , M 5 , 1 , 3 and 4 theoretically.To summarize the dynamics of system (2), we achieve the following theorem.Assume that the assumptions (A1)-(A3) hold.Then system (2) has five boundary equilibria: E 0 , E 1 , E 2 , E 3 , E 4 , and a unique positive equilibrium E * , in which E 0 , E 1 and E 2 are unstable, E 3 , E 4 are asymptotically stable.Assume that the parameters in system (2) are defined the same as in Example 3.3.Then the conditions for Theorem 3.4 are satisfied, and there exists six equilibria: E 0 , E 1 , E 2 , E 3 , E 4 and E * .Moreover, equilibria E 0 , E 1 , E 2 and E * are unstable, equilibria E 3 and E 4 are asymptotically stable.Furthermore, we can see that the assumptions (A1)-(A3) hold, and the equilibria E 0 , E 1 , E 2 are unstable, while the equilibria E 3 and E 4 are asymptotically stable, which are shown in Figure4.The theoretical results in Theorem 3.4 are thus numerically verified.
* can not be determined theoretically in this paper, since we only know the sign of M 1 .Nevertheless, the signs of M 2 , M 3 , M 4 , M 5 and 1 , 3 , 4 are definite under a specific parameters' setting.In the following, we present a numerical example to demonstrate the theoretical results in Theorem 3.4 and explore the stability of E * in a given parameters' setting.Figure 4. Example 3.4: Let the parameters α, d 0 , d 1 , μ, c, a, δ, m be specified the same as in Example 3.3.Then we have κ(m) ≈ 0.21046512 and κ(m) ≈ 0.27395349.Select β = 0.25 ∈ (κ(m), κ(m)).Then M 5 ≈ −0.00028413475 < 0, suggesting that the characteristic equation has at least one root with a positive real part.Thus the equilibrium E * is unstable, which is consistent with Figure 4.