Steady State Stability Analysis of a Chagas Disease Model

Chagas disease is caused by the parasite Trypanosoma cruzi, which is spread primarily by domestic vectors in the reduviid family. This work presents a model of the dynamics of Chagas disease in a rural village. The model consists of a nonlinear delay logistic-type differential equation for the total population of vectors and three nonlinear differential equations for the populations of infected vectors, infected humans, and infected domestic mammals. Steady state solutions for the model are derived and analyzed. Stability numbers are provided along with conditions for local asymptotic stability and partial results for global asymptotic stability. Numerical simulation results are presented, verifying the theoretical results. Keywords-Chagas disease, epidemic dynamics, delay logistic model, steady states, nonlinear dynamical system, 92D30, 92D25, 93C15, 34D20.


I. Introduction
Chagas disease is caused by the parasite Trypanosoma cruzi.It leads to organ deformity and early death in one third of the 8-10 million individuals infected throughout Latin America [24], [30].T. cruzi is primarily spread by reduviid vectors with some species particularly specialized in domestic infection cycles.Information on the relevant biological processes can be found in [10].Additionally, the parasite can be transmitted congenitally, orally, and through blood transfusions and organ transplants [4], [19], [21], [27].
Current control measures for the spread of the disease include improving the quality of housing and health care and treatment of homes with insecticide spraying [31].Additional control measures are treatment for acute Chagas disease and congenital transmission cases [17].However, because of the lack of an effective vaccine and the toxicity and questionable efficacy of the current drug treatments during the chronic stages of the disease, controlling the transmission of Chagas disease remains largely based on controlling the vector population and on blood-bank screening [12], [25], [29], [31].The spraying uses pyrethroid insecticides and has proven to be effective in retarding the spread of the disease and, in some cases, nearly eliminating the domestic insects [11], [22].Two mathematical models [33], [34], each consisting of four nonlinear differential equations, have been recently used to study the dynamics of Chagas disease transmission.Each system of equations models the total number of vectors, infected vectors, infected humans, and infected domestic animals in a rural village setting.The disease is transmitted by the carrier insects living in the village houses.In each of the models, the total number of vectors V = V(t) is modeled by a delay differential equation.In [34], a delayed Nicholson's blowfly-type term is used for the growth of vectors, whereas in [33], a delay logistic term represents the growth rate of the vectors.The former allows control of the maximum growth rate, and the latter allows control over the vector carrying capacity in the village houses and a vector death rate due to overpopulation beyond the carrying capacity.The delay logistic equations used in [33] and in this work are similar to the houseflies model presented in [35].However the model here differs in that the logistic term vanishes when the population is above the carrying capacity.For a survey of single species dynamics governed by delay differential equations, see [28] More recently, a mathematical model and simulations for the spread of Chagas disease in rural villages were presented in [9].The model is an enhanced version of that in [33], because it includes additional terms involving effects of congenital transmission in humans and domestic mammals and oral transmission in domestic mammals.Here, we present a steady state stability analysis for a special version of the model in [9]: the case where the model has constant coefficients, rather than the more general case of seasonal coefficients.
In Section 2, we present the model.The steady state stability analysis is presented in Section 3. In Section 4, numerical simulation results are reported and discussed in relation to the theoretical results.In Section 5, the stability analysis is summarized and unresolved issues are presented.

II. The Model
In this section, we present a model for Chagas disease dynamics in rural villages.We model the total number of vectors in the village houses using a delay differential equation, and we model the infected vectors, infected humans, and infected domestic mammals using ordinary differential equations.The model that we present has constant coefficients.However, a more general version of the model coefficients can be found in [9], where simulations are shown with temporal coefficients that are taken to be yearly periodic to represent seasonal changes.
We begin with the equation for the total number of vectors V = V(t) in the village houses at time t days.The rate of change in the total vector population depends on the egg hatching rate d h , which is a product of terms involving the following: the fraction of the vector population that can lay eggs, the number of eggs laid by an adult female per bite, the successful hatching of eggs after τ > 0 days (gestation time), the total blood supply b supply (measured in the unit of human factors), and the biting rate per human factor b h .The total rate of change also depends on the natural death rate d m of triatomines and the mortality rate d k due to overpopulation beyond the carrying capacity K.
Vector consumption by domestic mammals (referred to as dogs) is included in the model since mammals can become infected by predation on vectors [4], [19], [21], [27].Following the approach in [19], we use a Holling type II functional response for the consumption rate F(V): with units, vectors per dog per day.Here, E is the maximum number of vectors consumed by one dog per day, V min is the number of vectors found in the deep cracks of the village homes, and A + V min is the vector population at which dogs consume at the rate of E/2 vectors per day.Since we assume that dogs do not eat vectors in the deep cracks, the total vector population available for consumption is (V − V min ) + .Here and below, f + and f − denote the positive and negative parts of a function, respectively, so that f = f + − f − .Next, for simplicity, we assume that the number of dogs D in a village is constant.Using a delay logistic-type growth term, the rate of change of the total vector population in the village houses is In the first term, d h V(t − τ) is the rate of vectors per day that hatch at time t from eggs laid at time t − τ.The expression (1 − V(t − τ)/K) + represents the fraction of the food supply that was available to the female vectors at time t − τ.This assumes that if the vector population at time τ days prior to the current time t was above the carrying capacity, then no eggs were laid and, indeed, when K < V(t − τ) the first term vanishes.The natural death rate coefficient is d m and the death rate caused by overpopulation is described by The last term is the vector consumption term, which contributes to the removal of vectors.
Next, we consider the growth rates of the infected vectors V i , infected humans N i , and infected dogs D i .We denote the total number of humans by N and domestic non-mammals by C, and assume they are constant for simplicity.Following [10], we use 'human factors' as the unit for the total blood supply in the following way: each human represents one human factor, each domestic mammal represents d f human factors, and each domestic non-mammal represents c f human factors; then, the total blood supply is given by b supply = N + d f D + c f C. We denote by P NV and P DV the probabilities of a vector becoming infected by biting an infected human and an infected dog, respectively.Also, the number of bites per vector per day is given by b h • b supply , where b h is the same biting term used to define d h .Since the proportions of the bites that occur on infected humans and infected dogs are N i /b supply and d f D i /b supply , respectively, the growth rate of infected vectors is We assume that the natural death rate of the infected vectors is also d m , i.e., carrying the parasites does not affect their life span.Similarly, d k is the death rate due to lack of resources and space when the population size is above the carrying capacity.
The death rate of the infected vectors due to predation by dogs is G(V)DV i , where Collecting the terms above, we get the following rate equation for the infected vectors: We turn now to the rate equation for infected humans.When bitten by an infected vector, a susceptible human becomes infected with probability P V N .As before, each vector is biting at a rate of b h • b supply bites per day, and (N − N i )/b supply is the fraction of bites that are on susceptible humans.Thus, the growth rate of infected humans is given by b h P V N (N − N i ) V i .We assume that human reproduction is independent of infection status, since infected humans typically live beyond reproductive ages.Therefore, the assumption that N is constant implies that the rate of congenital transmission in humans is Here, T N i is the ratio, amongst babies born to infected mothers, of infected babies to the total number of babies.The mortality rate coefficients for infected humans, infected dogs, susceptible humans, and susceptible dogs are denoted by γ N i , γ D i , γ N s , and γ D s , respectively.
Therefore, the rate of change of infected humans is Under the assumption that both infected and uninfected dogs reproduce at the same rate, the rate of change of infected dogs is similar to that of infected humans, except that we must account for infection caused by the consumption of infected vectors.To this end, let P V D b be the probability that an uninfected dog becomes infected when bitten by an infected vector, let P V D c be the probability that an uninfected dog becomes infected after eating an infected vector, and let D − D i be the susceptible dog population.Then, the equation for the rate of change of infected dogs is where T D i is the probability that an infected dog passes the infection to its offspring congenitally, and the congenital transmission term is similar to the term in the infected humans equation because D is constant.
We note that domestic non-mammals, such as chickens, cannot become infected and thus are not considered as an infected population.However, they are available to the vectors as a blood meal and are included in the blood supply.
To complete the model, we assume that V min < K and prescribe the initial values of the respective populations: To summarize, we have the following system of equations:

III. Steady States and Stability Analysis
A steady state analysis of the model system (II.1)-(II.4) is presented in this section.Note that the vector equation is independent of the infected populations.Therefore, equation (II.1) is decoupled from the other equations in the system and can be analyzed separately.

A. Steady States and Stability of the Vector Equation
The steady states V of equation (II.1) are found by solving Clearly, V = 0 is a solution to equation (III.1),since when V = 0, (V − V min ) + = 0 and all other terms vanish, as well.We note that each steady state solution satisfies , and the sum of the other two terms on the right-hand side of (III.1) is negative.
and note that when R > 1, the vector growth rate exceeds the vector mortality rate.We next consider the case when 0 < V ≤ V min , so that (V −V min ) + = 0. Thus, equation (III.1)reads The positive solution is For R in this range, this is the only possible solution as we show below.
We now consider V min < V < K and then equation (III.1)can be written as (III.2) First, the factor (K R −V) in (III.2) must be positive for a solution V to be in (V min , K), which implies that In particular, this shows that there are no solutions such that Moreover, it is straightforward to see that the cubic polynomial in V on the right-hand side of equation (III.2) intersects the line (with positive slope KED/d h ) on the left-hand side of the equation at a unique positive value V in the interval (V min , K R ).Therefore, there is a unique steady state value V for R > K/(K − V min ).
We now turn to the stability of the steady states.We will linearize about each equilibrium and perform the stability analysis on the resulting equation, since it is well known that this is sufficient for studying the local stability of steady states for delay differential equations of this type [2], [3].
First we consider the V = 0 equilibrium.To this end, we linearize (II.1) about V = 0 to get Then, the local stability of the vector-free equilibrium follows from standard theory with V = 0 unstable for R > 1 and locally asymptotically stable for 0 < R < 1, see Theorem 4.7 from [32].
We next show that the vector-free equilibrium is, in fact, globally attracting for 0 < R < 1.To that end, first note that V(t) > 0 for all t > 0 when V(0) > 0. To see this, suppose V(t) < V min , which must occur if V is to reach 0. Then equation (II.1) gives dV/dt ≥ −d m V and V is bounded below by a decaying exponential function.Thus, V remains positive.Now, considering equation (II.1) again, we see that Considering (III.3) with equality instead of inequality, we see that 0 is globally asymptotically stable for 0 < R < 1, see [18], page 108.Thus, when V(0) > 0, inequality (III.3)implies that V converges to 0 as well.
Next, we turn to case with 1 < R < K/(K−V min ), where we have the positive equilibrium Then, so that, after simplifying, equation (II.1) becomes Next, let u(x) = y(x) − 1, so that the equation in terms of u is which can be linearized about u = 0 so that As before, the stability of the u = 0 equilibrium follows from standard results [32].Thus, when 1 < R < K/(K − V min ), the V = K R equilibrium is: locally asymptotically stable for 1 < R ≤ 3; locally asymptotically stable for R > 3 and d h τ < τ * h ; unstable for R > 3 and d h τ > τ * h , where .
Finally, we consider the case where R > K K−V min > 1 and V is the unique positive steady state solution that solves the cubic equation (III.2).We again let x = d h t and x = d h τ, and define Then, equation (II.1), in terms of y 1 , becomes And, the solution V ∈ (V min , K R ) to (III.1) is normalized to the solution y 1 = 1 in (III.4).
Next, to linearize at zero, we let u 1 (x) = y 1 (x) − 1, so that equation (III.4) in terms of u 1 becomes After some manipulation and the use of (III.2), this can be linearized to be where (III.5) Then the stability of the u 1 = 0 equilibrium follows as before [32].The stability of the positive steady state V ∈ (V min , K R ) is summarized below along with the other results from this section.
Theorem 1.The following hold for the vector equation (II.1).
1.The vector-free steady state, V = 0, exists for R ∈ (0, ∞), is globally asymptotically stable for 0 < R < 1, and unstable for R > 1. 2. When 1 < R < K/(K − V min ), there exists a unique positive steady state V = K R that is a. locally asymptotically stable for 1 < R ≤ 3, b.locally asymptotically stable for R > 3 and where

B. Steady State Analysis of Infected Populations
We now study the steady states of the system (II.2)-(II.4).For convenience, we use the following notation, We note that Γ N , Γ D > 0, since We first consider the case of R < 1, where V = 0 is the only vector equilibrium, and the system has only one biologically feasible steady state: the disease-free state (V i , N i , D i ) = (0, 0, 0).We note that there is an equilibrium (0, N i , D i ) with N i , D i > 0. However, this equilibrium violates N i < N, D i < D, and is therefore not biologically relevant.The disease-free state is globally asymptotically stable.Indeed, the Jacobian matrix at (0, 0, 0) is .
The eigenvalues at this steady state are λ 1 = −d m , λ 2 = −Γ N , and λ 3 = −Γ D , which are all negative, guaranteeing local asymptotic stability [1].The global convergence from positive initial conditions to (0, 0, 0) follows from Theorem 1.1 and inspection of the system (II.2)-(II.4).In particular, since V → 0, we see from (II.2) that V i → 0 as well.Then since V i → 0, (II.3) and (II.4) Now, we turn to the case when R > 1, where there are two steady state solutions for V : V = 0 and V > 0, where In what follows, we assume that V is positive, because the vectorfree equilibrium is unstable for R > 1.
For the sake of simplicity, we redefine The steady states of (II.2)-(II.4) are the solutions to the following system, where V is the solution of (III.1): (III.8) The system (II.2)-(II.4)clearly has a disease-free steady state, and we will later show that it has another positive steady state: (i) The disease-free state: Let us first consider the stability of the diseasefree steady state.The Jacobian matrix for the system (II.2)-(II.4) at the disease-free steady state is , where The characteristic equation for the eigenvalues λ of A can be simplified to the following: where It is reasonable to assume that d m > Γ D > Γ N > 0, since the half-life of vectors is smaller than that of dogs, which is smaller than that of humans.Also, a 0 12 , a 0 21 , a 0 13 , a 0 31 > 0, so that Γ N < B Γ /B < Γ D .From this, we see that the cubic polynomial C(λ) on the left-hand side of equation (III.9) has three distinct roots, − d m < −Γ D < −Γ N , all of which are negative, and the line L(λ) on the right-hand side of (III.9) has one root, −B Γ /B, in the interval (−Γ D , −Γ N ).Therefore, equation (III.9) has three real solutions with at least two being negative.The other solution will be negative if and only if L(0) < C(0), which is equivalent to With this, we define a stability number for the system of infected individuals as We conclude that the disease-free equilibrium is locally asymptotically stable when Π 0 < 1 and unstable when Π 0 > 1.
We now show that there is a positive endemic steady state corresponding to V > 0. It follows from (III.7) and (III.8) that For x ≥ 0, note that the expressions under both radicals above are nonnegative if which is equivalent to , (III.12) and which is equivalent to Inequalities (III.12) and (III.14) are reasonable for parameter sets found in the literature.In particular, the right-hand sides of both inequalities are greater than or equal to 0.5, while the data in [21] indicates that T N i is between 0.02 and 0.1.Therefore, we assume that the two inequalities (III.12) and (III.14)hold.Under these assumptions it also follows that y ± , z ± > 0.
Next, since y = N i ,we must have that y ± ≤ N, which is equivalent to Since y + violates (III.11), the only possible solution is y − which does, in fact, lie in (0, N] under the parameter assumptions and inequality (III.12).Similarly, z + ≤ D violates (III.13), and z − ∈ (0, D] under the parameter assumptions and inequality (III.14).For notational convenience, we drop the superscripts on the positive solutions y − and z − .Then, (III.6) can be written as where .
Since y, z > 0, it follows that F is a positive function on [0, ∞).Furthermore, inequalities (III.11) and (III.13)imply that F is an increasing function on [0, ∞).Therefore, xF (x) is an increasing function on this interval.Moreover, xF (x) > xb h on [0, ∞).Thus, there exists a unique solution x to (III.15) in the interval (0, V), and there is exactly one positive endemic steady state: (x, y, z) in (0, V) × (0, N] × (0, D]. We now analyze the stability of the unique positive equilibrium point, (x, y, z), with V > 0. The associated Jacobian matrix of (II.2)-(II.4) has column vectors For convenience in what follows, we denote the entries of the Jacobian matrix above by (a i j ), where i, j = 1, 2, 3. Notice that the parameter assumptions and inequalities (III.11) and (III.13)imply that all three diagonal entries a ii , i = 1, 2, 3, are negative, while all other nonzero entries are positive.The characteristic equation for the eigenvalues λ of the Jacobian matrix can be simplified to Then, B, B Γ > 0 and we have the following bound: Therefore, the line 3 (λ) on the right-hand side of equation (III.16) has a negative λ-intercept that lies between a 22 and a 33 , which are two roots of the cubic polynomial c(λ) on the left-hand side of the equation.This implies that all three roots of the characteristic polynomial are real numbers and that at least two of them are negative.The third root is negative if and only if 3 (0) < c(0), which is equivalent to B Γ < −a 11 a 22 a 33 .
We now define another stability number and note that when conditions (III.11) and (III.13)hold, the endemic steady state is locally asymptotically stable for Π < 1 and unstable for Π > 1.We point out that Π depends on the equilbria values (x, y, z) and that Π 0 is, in fact, Π evaluated at (0, 0, 0).However, inequalities (III.11) and (III.13) are not required conditions in the disease-free stability analysis.Finally, we note that the stability numbers depend on all the parameters of the system.The results of this section are summarized below.

IV. Numerical Results
In this section, we verify the results in Theorems 1 and 2 using a numerical algorithm that was coded in Mathematica [36].Baseline parameter values, which are our best estimates, can be found in Table I.In order to satisfy the hypotheses of the theorems, we tweaked parameter values as described below.In each case, where certain parameters are changed, all other parameters are kept at their baseline values.
A. Theorems 1.1 and 2.1 In order to obtain R < 1, we choose d h = 0.95d m .Fig. 1 depicts the global asymptotic stability of the vector-free and disease-free equilibria in this case.

B. Theorem 1.2
Since K/(K − V min ) ≈ 1.04 with the baseline parameters, we choose d h = 1.03d m to satisfy the hypotheses of Theorem 1.2a.Fig. 2 depicts the local asymptotic stability of the equilibrium V = K R .
To satisfy the hypotheses of Theorem 1.2b, we set V min = 5000 and d h = 3.05d m .The equilibrium V = K R is again locally asymptotically stable in this case, as verified by Fig. 3.
For Theorem 1.2c, we decrease the carrying capacity K to 100 vectors per house while keeping V min = 5000 and d h = 3.05d m .The oscillatory solution and instability of the equilibrium V = K R are shown in Fig. 4. The vector population under the hypotheses of Theorem 1.2a with a large and small initial population.The vector population under the hypotheses of Theorem 1.2b with a large and small initial population, shown over two and fifty years.

C. Theorem 1.3
Using the baseline parameters, we have R > K/(K − V min ) and A ≈ −L.Slightly lowering the carrying capacity to 400 vectors per home satisfies the hypotheses of Theorem 1.3b, and the local asympotic stability of the positive equlibrium V ∈ (V min , K R ) can be seen in Fig. 5 for this case.The baseline parameters without any changes fulfill the hypotheses of Theorem 1.3c and, in this case, the local asympotic stability of the positive equlibrium can be seen in Fig. 6.Despite many attempts, we were unable to change the parameters in a way that satisfied the hypotheses of Theorems 1.3a and 1.3d.It seems these cases are unlikely to

D. Theorems 2.2 and 2.3
Using the baseline parameters, we have R > 1, Π 0 > 1, and Π < 1.From Theorems 2.2 and 2.3 we expect an unstable disease-free equilibrium and a locally asympototically stable endemic equilibrium.This is verified in Fig. 6.

V. Conclusion
This work presents a local and partial global stability analysis for the steady states of a model for the dynamics of Chagas disease.The model consists of a delay logistic-type equation for the total number of vectors in village houses and nonlinear differential equations the populations of infected vectors, infected humans, and infected domestic mammals.In Section 2, a model for Chagas disease dynamics within village houses is presented.
In Section 3, the steady states are analyzed.The vector equation is decoupled from the other three equations and thus analyzed separately.The birthdeath ratio, R, controls the stability of the vector equilibria with the vector-free state being globally asymptotically stable for R < 1 and unstable for R > 1, as expected.
When R > 1, there is also a positive vector equilibrium whose stability is classified in Theorem 1.The conditions for local asymptotic stability are complex and depend on the parameters of the system.When Moreover, when R > K/(K − V min ), there is a unique positive equilibrium between V min and K R .In practice, this equilibrium was found to be asymptotically stable across a broad range of reasonable parameter sets.In fact, we were not able to find a parameter set that would produce an unstable positive equilibrium when R > K/(K − V min ).
In Section 3.2, we analyzed the stability of the steady states of the infected populations and report the results in Theorem 2. When R < 1, the diseasefree state is the only equilibrium and it is globally asymptotically stable.When R > 1, there are two equilibria, the disease-free state and endemic state.The stability of these equilibria is controlled by the stability numbers Π 0 and Π, respectively, where each depends on the parameters of the model.In particular, the disease-free state is locally asymptotically stable for Π 0 < 1 (unstable for Π 0 > 1) and the endemic state is locally asymptotically stable for Π < 1 (unstable for Π > 1).It is challenging to determine the biological meaning of these stability numbers, however it appears that for realistic parameter sets, Π 0 > 1 and Π < 1.That is, for the simulations we performed across a wide variety of reasonable parameter sets, the endemic state was always locally asymptotically stable.These results are consistent with the endemic nature of Chagas disease in Central and South America.The relationship between Π 0 and Π is of interest and open to further study.We performed many more numerical simulations than are presented here and the results suggest that our local stability conditions may be global ones.In the future, we intend to explore this possibility.
Fig.2.The vector population under the hypotheses of Theorem 1.2a with a large and small initial population.
Fig.3.The vector population under the hypotheses of Theorem 1.2b with a large and small initial population, shown over two and fifty years.

Fig. 4 .
Fig.4.The vector population under the hypotheses of Theorem 1.2c shown over three and fifty years.

Fig. 5 .
Fig.5.The four populations under the hypotheses of Theorem 1.3b.

Fig. 6 .
Fig.6.The four populations under baseline parameters with large and small initial populations.
unstable for A < −L and d h τ > τ * * h , where A and L are from (III.5) and τ * * h there exists a unique positive steady state V ∈ (V min , K R ) that is a. unstable for A > L, b. locally asymptotically stable for −L < A < L, c. locally asymptotically stable for A < −L and d h τ < τ * * h , d.

TABLE I The
model parameters and the baseline simulation values