TRANSMISSION DYNAMICS OF A CHAGAS DISEASE MODEL WITH STANDARD INCIDENCE INFECTION

In this paper, an insect-parasite-host model with Ricker’s type re-production of triatomines and the standard incidence rate of the interaction between insects and hosts is formulated to study the transmission dynamics of Chagas disease. Two thresholds of the ecological basic reproduction number of triatomines and the epidemiological basic reproduction number of Chagas disease are derived, which determine the dynamics of this model. As a result, the existence of equilibria and the local/global stabilities of the equilibrium are accordingly obtained. Moreover, backward bifurcation, forward bifurcation and saddle-node bifurcation are also shown analytically and numerically. Biologically speaking, Chagas disease may undergo outbreak if the number of bites of per triatomine bug per unit time or the transmission probability from infected bugs to susceptible competent hosts per bite increase.


Introduction
Chagas disease, known as American trypanosomiasis, which is spread through bites by triatomine bugs infected with T. cruzi or transmission from mother to child.It is early discovered in 1908 by Doctor Chagas and is considered as a neglected disease while still attracts public health attention.It is reported that about 13% of the Latin American population is at risk of Chagas disease infection [23], and it also has been observed considerable spread of this disease due to the convenient transportation and the widespread development of globalization [6,25,[27][28][29].The disease is transmitted by a protozoan parasite called by Trypanosoma cruzi (T.cruzi) by invading the lymphatic system and the blood stream.It is mainly prevalent in Central and South America, such as Argentina, Bolivia, Brazil, Chile, etc.The disease can cause serious clinical symptoms involving cardiomyopathy and heart disease, which may cause many infected people to die due to these symptoms [4].People infecting with T. cruzi generally passes through acute and chronic stages.
The former includes fever, facial edema, anemia and so on and the latter often has myocarditis, cerebral embolism, sudden death and so on [21].There are an estimated 6-8 million persons in the Americas with Chagas disease, causing a burden of 29,000,000 DALYs and a health care cost of 24.73 billion.
Recently, Wu et al. [36] formulated a new model (1) considering the Ricker's type function growth of triatomine bugs and pathogenic effect on triatomine bugs to study the dynamics of interaction between triatomine bugs and T. rangeli, which is relevant to Chagas disease.The model is S v (t) = r(S v (t) + θI v (t))e −σ(Sv(t)+Iv(t)) − βv S v (t)I h (t) − β c S v (t)I v (t) − µ 2 S v (t), where the population is divided into four components: susceptible and infected competent hosts, susceptible and infected triatomine bugs, denoted by S h , I h , S v , I v , respectively.Λ h is the constant recruitment rate of susceptible competent host per unit time.The transmission rate from infected bugs to susceptible competent hosts is denoted by βh = ba Nc+αNq , where b is the transmission probability from infected bugs to susceptible competent hosts per bite, a is the number of bites of per triatomine bug per unit time, α is the biting preference of quasi-competent hosts to competent hosts, N c is the total number of competent hosts, N q is the total number of quasi-competent hosts.The transmission rate from infected competent hosts to susceptible bugs is denoted by βv = ca Nc+αNq , where c is the transmission probability from infected hosts to susceptible triatomine bug per bite.The total infection rate through co-feeding transmission between susceptible and infected bugs is β c , which is transmitted by both the competent and quasi-competent hosts.The Ricker's type function b(x) = rxe −σx was chosen to model the reproduction rate of R. prolixus.Integrating the pathogenic effect, the growth rate of triatomine bugs is modeled as r(S v + θI v )e −σ(Sv+Iv) , where r is the maximal number of offsprings that a triatomine bug can produce per unit time, and θ ∈ [0, 1] is the reproduction reduction of bugs due to the pathogenic effect of T. rangeli on bugs, σ is the densitydependency strength measuring the reproduction of bugs.µ 1 and µ 2 are the natural death rates of competent hosts and triatomine bugs, respectively.d v is the death rate of infected vectors induced by pathogenic effect.
Model (1.1) was used to develop the systemic and co-feeding transmission routes among vectors and hosts, and two thresholds was derived to characterize the dynamical behavior of this model.Interestingly, sustained oscillations of model solutions were observed numerically by altering parameter d v and θ and the results have shown that the oscillation amplitude will be larger if d v is larger or θ is smaller.
As we know, Trypanosoma rangeli (T.rangeli) is a kind of parasite which is pathogenic to some vector species including triatomine bug, although it can infect The rest of the paper is organized as follows.In the next section, a new model with the death rate of the infected hosts, the standard incidence rate are formulated.The corresponding basic reproduction number is given.In section 3, the existence, and local and global stability of equilibria are obtained.In section 4, the forward bifurcation and the backward bifurcation are derived by using the center-manifold theorem.Numerical simulations are also performed to illustrate the obtained results.Finally, conclusion and discussion are also given.

Model description
Motivated by Wu et al. [36], in this paper, a new chagas model involved by the effect of the T. cruzi-induced death rate of hosts and the standard incidence rate between hosts and triatomine bugs is formulated.We assume the generation rate of triatomine bugs follows Ricker's type function instead of the logistic growth, and further study the dynamical behavior of triatomine-host transmission.The interaction among hosts, T. cruzi parasite and triatomine bugs are described by the following model with Ricker's type reproduction function of triatomine bugs: where all the parameters are non-negative, and its biological meanings and ranges are given in Table 1.
Let N h (t) = S h (t) + I h (t).Summing up the first and the second equations of system (2.1), we obtain N h (t) = Λ h − µ 1 N h − dI h , then we have the feasible region of system (2.1) is µ1 .Summing up the third and the fourth equations of system (2.1) and using so that we can immediately obtain lim sup which indicates the boundedness of S v and I v mentioned above.Now, let us calculate the basic reproduction number of system (2.1) by the method of Dreessche and Watmough [8]. Define Then the T. cruzi basic reproduction number of system (2.1) is given by the spectral radius of the next generation matrix F V −1 , which is where R v = r µ2 .

Existence and stability of equilibria
3.1.Existence of equilibria Obviously, system (2.1) has a vector-free equilibrium E 0 (S 0 h , 0, 0, 0) and a parasitefree equilibrium E S (S 0 h , 0, S 0 v , 0).Now we turn to analyse the existence of positive equilibria

1), its coordinates satisfy
, and it is easy to see that I * h , S * v and I * v are greater than zero if S * h is greater than zero.The coordinate of S * h should be the positive root of the following cubic equation: where For simplicity, we assume that A = 0, then we can rewrite equation (2.1) as where .
The derivative of equation (3.4) is which has two real roots . According to system (2.1), we have and more details are provided in Appendix I.
Then the following theorems are obtained by applying Cardano formula [14,37]: 1) has at most two positive equilibria.Moreover, system (2.1) has (1) two positive equilibria if (2) one positive equilibrium if (3) no positive equilibrium in other cases.
Proof.Obviously, the number of positive equilibrium of system (2.1) is corresponding to the number of different positive roots of cubic equation (3.4).We can see that f (0) = B 3 > 0 and f (±∞) = ±∞ when acµ 1 > dµ 2 , which implies that system (2.1) has at most two positive equilibria.
When condition (1) holds, namely, ∆ 3 < 0, B 1 ≤ 0, B 3 > 0 or ∆ 3 < 0, B 1 > 0, B 2 < 0, B 3 > 0, we know that the cubic it is impossible for equation (3.4) to have two negative roots by Veda's theorem, which means that there is at least one positive root η + of equation (3.4).By case III3 and III4 in Table 2 in Appendix I and B 3 > 0, we conclude that f (S h ) has two different positive roots if ∆ 3 < 0, which is described in Figure 1 (a) Then we proved case (1).
By cases II2 and II5 in Table 2 in Appendix I, we can know that it is a root of multiplicity 2 when equation (3.4) has one positive root and B 3 > 0. It is described in Figure 1 (b).This completes the proof.Theorem 3.2.If acµ 1 < dµ 2 , system (2.1) has at least one positive equilibrium and at most three positive equilibria.Moreover, system (2.1) has (1) three positive equilibria if Proof.The number of positive equilibrium of system (2.1) is corresponding to the number of positive root of cubic equation (3.4).When acµ 1 < dµ 2 , f (0) = B 3 < 0 and f (±∞) = ±∞, f (S h ) has at least one positive equilibrium and at most three positive equilibria.When ∆ 3 < 0, B 1 < 0, B 2 > 0, B 3 < 0, there are three different positive roots of the cubic equation (3.4) by case III4 of Table 2.When condition (1) holds, f (S h ) has three different positive roots by Descartes' rule of sign which is described in Figure 2 (a).
If f (S h ) has two different positive roots, one of the two roots is of multiplicity 2, which implies that f (η + ) = 0 or f (η − ) = 0, see Figure 2 (b).By case II4 of Table 2, we know that f (S h ) has two different positive roots if ∆ 3 = 0, B 1 B 2 < B 3 < 0. In addition, we add a constraint which is B 2 1 − 3B 2 > 0 when f (S h ) has a root of multiply 3 to ensure that f (S h ) has two different positive roots.Then case (2) is proved.From previous analysis, we know that f (S h ) has at least one positive equilibrium due to B 3 < 0, which means that there is one positive equilibrium in other cases in Table 2 when B 3 < 0. Thus, we complete the proof.(2) one positive equilibrium if ab ln R v < Λ h σ; (3) no positive equilibrium in other cases.
Proof.If acµ 1 = dµ 2 , then equation (3.2) becomes the following quadratic equation: where It is easy to see that D and ∆ are greater than zero.
If ab ln R v > Λ h σ, then B > 0. Therefore, there are two positive equilibria as long as C < 0, otherwise, there is no positive equilibrium.We complete the proof of case (1).
If ab ln R v < Λ h σ, then B < 0. According to the diagram of the quadratic equation (3.9) under these conditions, we can know that there is only one equilibrium.We complete the proof of case (2).
Theorem 3.5.The vector-free equilibrium E 0 (S 0 h , 0, 0, 0 Proof.Let N v = S v + I v and sum up the third and fourth equations of system (2.1), we have In view of N v (0) = S v (0) + I v (0), we have It means that the solutions of S v and I v with feasible initial condition tend to zero if R v < 1.For N h , it is cooperative with its positive invariance set [0, Λ h µ1 ].We also know that E 0 is a unique equilibrium of system (2.1) when R v < 1.Thus, we know that E 0 is globally asymptotically stable if R v < 1, and unstable if R v > 1.
Proof.System (2.1) has a parasite-free equilibrium E S (S 0 h , 0, S 0 v , 0) when R v > 1.It is easy to obtain the following Jacobian matrix J(E S ) of system (2.1) at E S , Then we have the following characteristic polynomial where The eigenvalues of Jacobian matrix (3.12) are denoted by where ., E S is locally asymptotically stable.
If R 0 = 1, we have λ 3 = 0, so that we need to illustrate the type of E S .Firstly, let to move the equilibrium E S to the origin, then system (2.1) becomes By using the transformation where Then we obtain which shows that E S is a saddle-node point when R 0 = 1.

Stability of positive equilibrium
In this subsection, we will study the local stability of positive equilibrium E * (S * h , I * h , S * v , I * v ).The Jacobian Matrix of system (2.1) at E * is It is easy to obtain the characteristic polynomial of J(E * ): ) where Then by using Routh-Hurwitz criteria, we have the following theorem: (1) acµ 1 > dµ 2 , and (2) acµ 1 = dµ 2 , and ab ln R v < Λ h σ.

Backward bifurcation and forward bifurcation
Now, we will study the forward bifurcation and the backward bifurcation of system (2.1) by calculating the direction of transcritical bifurcation point in the following. .
Proof.The Jacobian matrix at E S (S 0 h , 0, S 0 v , 0) of system (2.1) is Then the corresponding eigenvalues of the Jacobian matrix J(E S ) are, respectively, , Because of the biological significance, we know that all parameter values are nonnegative.Obviously, eigenvalues λ 1 and λ 2 are negative.
Let R 0 = 1, then we have Replacing c in λ 3 and λ 4 with c * , we obtain that λ 3 = 0, λ 4 < 0. From Theorem 3.6, we know that parasite-free equilibrium E S (S 0 h , 0, S 0 v , 0) of system (2.1) is unstable if c > c * , and locally asymptotically stable if c < c * which means that c = c * is a bifurcation value.
Through calculation, we get a right eigenvector u and a left eigenvector v associated with the zero eigenvalue as follows: By using orthogonal condition < u, v >= 1, one yields .
Next, set then system (2.1) has the form The bifurcation coefficients in system (2.1) at E S are as follows: Obviously, since v 1 = v 3 = 0, we only consider the cross derivatives of f 2 and f 4 in system (4.3) at E S .Then we can obtain non-zero terms as follows: The corresponding values of a * and b * are as follows: From [37,38], we realized that the signs of a * and b * determine the local dynamical behaviour around parasite-free equilibrium E S of system (2.1).It is easy to see that b * > 0 since S * h are great than zero.There are an unstable equilibrium exhibiting backward bifurcation near E S if a * > 0, and a locally asymptotically stable equilibrium undergoing forward bifurcation if a * < 0. Then we will discuss the sign of a * to study the dynamics of system (2.1).Define Then we can draw conclusions that system (2.1) undergoes a backward bifurcation if Λ h < Λ * h , i.e., a * > 0; If Λ h > Λ * h , i.e., a * < 0, system (2.1) undergoes a forward bifurcation.This completes the proof.

Numerical simulations
In this section, we will illustrate numerically the results of system (2.1) with the following parameter values in Wu et al. [36] by using Matlab and Auto07P [7]:   Using the above parameter values, it is easy to calculate that there is a vectorfree equilibrium (2000, 0, 0, 0), a parasite-free equilibrium (2000, 0, 5971.437,0) and a parasite-positive equilibrium (16.081, 661.306, 167.827, 5803.610).Firstly, we simulate the backward and forward bifurcation to verify the results in the previous sections.We choose transmission probability from infected hosts to susceptible triatomine bug per bite (c) as the primary bifurcation parameter and keep the other parameters fixed as (5.1), then we get the one-parameter bifurcation diagram, shown in Figure 3.In Figure 3 (a), there is a transcritical bifurcation point T C(4 × 10 4 , 0, 5.97144 × 10 3 , 0) showing backward bifurcation with Λ h = 100 when c = 0.0193049, and in Figure 3 (b), there is a transcritical bifurcation point T C(2 × 10 3 , 0, 5.97144 × 10 3 , 0) showing forward bifurcation with Λ h = 5 when c = 0.0009652.A numerical solution of system (2.1) tending to a stable equilibrium when the time tends to infinity with initial solution (S h , I h , S v , I v ) = (172, 609, 21, 576) is shown in Figure 4, where Λ h = 5, r = 0.0274, σ = 0.0002, µ 1 = 0.0025, µ 2 = 0.0083, d = 0.0083.Secondly, when the number of bites of per triatomine bug per unit time (a) is used as the primary bifurcation parameter, then we obtain the bifurcation diagram, shown in Figure 5.There is a transcritical bifurcation point T C(2000, 0, 5971.4377253,0) showing forward bifurcation with a = 0.02663.It is obvious that infected hosts (I h ) and infected vectors (I v ) will increase as the parameter a increases, which indicates that chagas disease may go through outbreak.
Thirdly, the transmission probability from infected bugs to susceptible competent hosts per bite (b) is used as the primary bifurcation parameter, we obtain the bifurcation diagram, shown in Figure 6.A transcritical bifurcation point T C(2000, 0, 5971.4377253,0) shows forward bifurcation as b = 0.0001182.The infected hosts (I h ) and infected vectors (I v ) will increase as the parameter b increases which indicates that chagas disease may go through outbreak.
Finally, we take the T. cruzi-induced death rate of hosts (d) as the primary bifurcation parameter, the bifurcation diagram shown in Figure 7 is obtained.The number of infected hosts (I h ) and infected vectors (I v ) will decrease as the parameter Biologically speaking, chagas disease may undergo outbreak if the number of bites of per triatomine bug per unit time (a), the transmission probability from infected bugs to susceptible competent hosts per bite (b) increase while chagas disease may disappear if the T. cruzi-induced death rate of hosts (d) decreases.

Discussion and conclusion
In this paper, a chagas model (2.1) with Ricker's type function and the standard incidence rate is investigated by using the dynamical system approach, and the existence and stability of equilibria are obtained, the related bifurcations are given, which tells that the positive equilibrium of model (2.1) is locally asymptotically stable when the number of bites of per triatomine bug per unit time (a), the transmission probability from infected bugs to susceptible competent hosts per bite (b) and T. cruzi-induced death rate of hosts (d) increase.Once the standard incidence rate is considered, the positive equilibrium of system (2.1) will be stable, however, it still tells us that the outbreak of infected hosts infected vectors without sustained oscillations.

Figure 1 .
Figure 1.The number of root for f (S h ) = 0 when acµ1 > dµ2.(a) Two different positive roots.(b) One positive root.

Figure 3 .
Figure 3. (a) Backward bifurcation diagram of system (2.1) showing c vs S h with Λ h = 100; (b) Forward bifurcation diagram of system (2.1) showing c vs S h with Λ h = 5.

Figure 5 .Figure 6 .
Figure 5. One-parameter bifurcation diagram of system (2.1) with respect to a. (a) a vs.I h ; (b) a vs. Iv.

Table 1 .
Parameter description of system(2.1) T ) for some T > 0. It remains to prove the boundedness.We know that S h +I h ≤