Analysis of a nonlinear system for community intervention inmosquito control

Non-linear difference equation models are employed in biology to 
describe the dynamics of certain populations and their interaction 
with the environment. In this paper we analyze a non-linear system 
describing community intervention in mosquito control through 
management of their habitats. The system takes the general form: 
 
 
 $x_{n+1}= a x_{n}h(p y_{n})+b h(q y_{n})$ 
n=0,1,... $y_{n+1}= c x_{n}+d y_{n}$ 
 
 
 
 
where the function $h\in C^{1}$ ( [ $0,\infty$) $\to $ [$0,1$] ) satisfying 
certain properties, will denote either $h(t)=h_{1}(t)=e^{-t}$ 
and/or $h(t)=h_{2}(t)=1/(1+t).$ We give conditions in terms of 
parameters for boundedness and stability. This enables us to explore 
the dynamics of prevalence/community-activity systems as affected by 
the range of parameters.

1. Introduction. Nonlinear difference equation systems are utilized in biology in order to describe the dynamics of certain populations. The seminal work of May et al. [16] offers eloquent examples in this regard. Studies of global stability, boundedness, and persistence in nonlinear discrete models of order greater than or equal to one have been carried out in a number of books and monographs: Agarwal [1], Kocic and Ladas [13], Caswell [4], Elaydi [8], Cushing [5] as well as articles: Kuang and Cushing [12], Allen et al. [2], Kulenovic and Nurkanovic [15], Freedman and So [9], Basson and Fogarty [3] etc. For one dimensional population models (first-order difference equations), some general global results have been obtained by Cull ([6], [7]), Rosenkrantz [17], and Huang [11]. Global stability theorems for two dimensional or higher-order systems are less abundant and a good survey of some partial results is Kulenovic and Merino [14] pp. 75-200.
In this paper we describe an exponential as well as a rational model for community intervention in mosquito control. We analyze the dependence between the level of consciousness of the general public and the dynamics of the breeding sites. The aim of our analysis is twofold. First we analyze the global stability behavior of the models in a more general context and give mathematical results that complement the studies mentioned in the introductory paragraph. In addition, we investigate numerically the regions of cyclic behavior. Second, we discuss how the models proposed here are relevant to community intervention in mosquito control.
The breeding sites of urban mosquitoes (for example, those that transmit the dengue virus) are formed when people pollute the environment with vases, old tires, cans, and other containers that are filled during the rainy season with water and where mosquitoes can lay their eggs. Breeding sites are often removed by public action, which requires not only awareness but also organization. The point of view we take in the present paper is that individual behavior determines the "birth rate" of breeding sites, while collective action predominates in the removal or "death rate." This leads to a mathematical model in the form of the system of difference equations    x n+1 = ax n h(py n ) + bh(qy n ) n = 0, 1, . . . , y n+1 = cx n + dy n (1) where x and y are dynamical variables quantifying the number of breeding sites and the "level of consciousness," respectively. As we shall see, the level of consciousness is an underlying variable which is defined through an assumed relationship between the fractions f C n = h(py n ) ∈ (0, 1] and f I n = h(qy n ) ∈ (0, 1]. The former quantity denotes the fraction of breeding sites that survive after the community intervention, whereas the latter is the fraction of breeding sites that are still newly created by individuals with a certain level of consciousness. The assumed relationship is that the fractions f C n and f I n depend in a similar fashion on a common dynamical variable y n , except perhaps for different sensitivities that can be tuned by modifying the parameters p and q. The number of sites newly created as a result of everyday activity by individuals with no social consciousness is denoted by b > 0. This number is expected not to vary on a weekly basis, save for cases of major technological or social advances that result in lifestyle changes. An increase in the level of individual consciousness (due, for example, to educational campaigns) produces a reduction in the number of newly created breeding sites to the value bf I n = bh(qy n ). Notice that the number of breeding sites that are newly created does not depend on the variable x n , which is the actual number of breeding sites. Through this assumption, we acknowledge that the individual possesses neither the means nor the desire to evaluate the magnitude of the infestation. The observation of breeding sites requires organization and community intervention. The latter can manifest itself indirectly, through educational campaigns, or directly, through removal of a portion of the breeding sites. The first effect has already been accounted for by assuming a reduction in the number of newly created sites. The second effect is a reduction of the actual number of breeding sites to the value ax n h(py n ), where the coefficient a ∈ (0, 1) denotes the fraction of the breeding sites that would naturally survive to the next week if the community intervention were absent. In the present model, the consciousness y n+1 increases linearly with the number of breeding sites x n from the previous week, through a constant c ≥ 0. The parameter c ensures a feed-back mechanism between the number of breeding sites and the level of consciousness. Of course, there is always a remainder of awareness from the previous week (this is a memory effect), which is quantified by dy n , with 0 < d < 1. Thus, the parameter (1−d) represents the rate of erosion of the level of consciousness. The linear dependence between y n+1 and x n is not a crucial limitation of the present model for, in fact, the relation between the variable y n and the fractions f C n and f I n is contingent upon the specific form of the function h. The analytical form of the function h is not exactly known. However, the arguments presented in the preceding paragraphs show that h must satisfy the following mathematical properties (P1): i) h ∈ C 1 ([0, ∞) → (0, 1]), ii) h is non-increasing, and iii) h(0) = 1 and lim y→∞ h(y) = 0. The function h is assumed to be nonincreasing because the role of consciousness through activity is to reduce the survival and creation of breeding sites. If y = 0 (no consciousness at all), we naturally have h(0) = 1, which means that the death of the breeding sites is attributed to natural causes alone, as modeled by the factor a. The property lim y→∞ h(y) = 0 claims that a sufficiently high level of consciousness leads to the death of all breeding sites within the current week. The variables and parameters entering system (1) are summarized in Table 1. Despite all our explanations, the reader may still feel unsure about the significance of the level of consciousness y n , which appears to be an implicitly defined variable connecting the fractions f C n and f I n . Let us remember that f C n is the fraction of breeding sites that survive after community intervention in week n, whereas f I n is the fraction of breeding sites that are newly created in spite of the individual's efforts to comply with the educational campaigns. In principle, these fractions can be determined through observation or experiment. Clearly, f C n ∈ (0, 1) is a measure (level) of community consciousness. Nevertheless, any monotonic function of f C n is a measure of community consciousness, as well. After all, we only need to establish a criterion that allows us to compare the levels of community consciousness from week to week or from site to site (that is, an order relation). Likewise, f I n ∈ (0, 1) is a measure of individual consciousness and so is any monotonic function of f I n . In particular, because the function h −1 is monotonic, y C n = 1 p h −1 (f C n ) and y I n = 1 q h −1 (f I n ) are measures of the community and individual consciousness, respectively. We assume that, when formulated this way, y C n = y I n , that is to say that f C n and f I n are in the relationship specified by Thus, the common value y n = y C n = y I n is a measure of both the community and individual awareness levels (we generally call it level of consciousness) and represents a convenient way to write the equations defining the dynamics. Of course, we could have expressed everything in terms of the fraction f C n (or, according to equation (2), in terms of f I n ) but doing so results in a rather complicate-looking system: Even though the preceding system adequately positions the model in its eco-biological context, for mathematical purposes, it is more convenient to work with the system (1). Also notice that, if f C n+1 is accepted as the natural measure of the community consciousness for the next week, then the dependence with x n expressed by the second of the two equations in the previous system is no longer linear. As already explained, the number x n+1 of breeding sites for the next week decreases when the awareness (or consciousness) y n is high. The effect of consciousness is either a reduction in the number of sites ax n (here, 0 < a < 1) that naturally survive to the next week or a reduction in the number of sites that are newly created. In both cases, the reduction is described by the same functional dependence h(t) but with different sensitivity parameters p ≥ 0 and q ≥ 0, with p + q > 0. The individual and the community are expected to exhibit different sensitivities to the same problem. If α = q/p ≤ 1, let us agree that the community consciousness is stronger than the individual one, since f C n ≤ f I n . Although both the local and global analysis of system (1) is performed in terms of the unknown function h, we also consider two specialized functional forms of h: a form with a fast decay, which is specified by the exponential dependence h(y) = h 1 (y) = e −y , and a form with a slow decay, which is specified by the rational dependence h(y) = h 2 (y) = 1/(1+y). By choosing the two functional forms for h(y), we span a large pallet of possibilities and conclusions regarding system (1). Both the exponential and the rational models capture the reality that an increase in consciousness reduces the survival and formation of new breeding sites. The models are analyzed for boundedness and persistence. We find a tight attracting interval that proves instrumental in proving global attractivity results. We also give a sufficient condition for global asymptotic stability in terms of parameters for the rational system.
In order to reduce the number of parameters in our equations, we distinguish between the cases when p and c are both strictly positive and the case when at least one of them is zero. After rescaling the equations in order to reduce the original six-parameter model to a four-parameter model, we shall analyze the model for mathematical properties such as boundedness, persistence, and stability, all of which have implications with respect to the kind of measures that are to be taken by health officials. In section 1.1, we prove the existence and uniqueness of the positive equilibrium and we obtain a sufficient condition for local stability in terms of system parameters. The model for which the number of breeding sites x n+1 for the next week decreases exponentially with the increase in awareness y n will be called the exponential model. The rational model which is introduced in section 1.1 has the property that the decrease in x n+1 is modeled by the factors 1/(1+py n ) and 1/(1 + qy n ), respectively. Sections 2 and 3 deal with attracting intervals and global stability results. We only summarily touch upon the local stability analysis. Instead, we refer the reader to the extensive literature existing on this better-understood aspect of the problem (Kulenovic and Merino [14] pp. 98-100, Agarwal [1], Elaydi [8], and Kocic and Ladas [13]).
The mathematical results in this paper are relevant to community intervention in the sense that they tell us where to concentrate our efforts (individual or community action). The goal of a successful public health policy is the design of interventions such as the development of educational programs that will reduce the number of mosquito breeding sites. The dynamics of this model reveals the fact that we need to concentrate our efforts on the community awareness as expressed through the removal of breeding sites. Last but not least the mathematical results are used to show that the prevalence of breeding sites exhibit some periodic behavior in some parameter regions, instead of converging to the desired equilibrium at the very low level. We do not have numbers for the parameters, however we know which actions can increase or decrease them. For example, increasing c, means observing breeding sites and reporting on breeding sites. An increase in d, indicates a higher survival of the previous awareness, which can be augmented by social campaigns. The parameter p characterize the conversion of consciousness into action. To increase p, means organization such as people going out and cleaning out the breeding sites. (1). In this section, we perform a stability analysis on the general model (1) linking mosquito breeding sites x n and community intervention y n , in the case when p > 0, c > 0. We keep this assumption on the parameters p and c throughout sections 1 and 2. The degenerate cases that happen when at least one of the parameters p or c are zero are treated separately in section 3. By the change of variables y n = z n /p and x n = r n /(pc), we get rid of two of the parameters. Denote β = pbc and α = q/p. Clearly β > 0 and α ≥ 0. After rescaling, the system of equations (1) takes the form:

Local Stability Analysis of System
   x n+1 = ax n h(y n ) + βh(αy n ) n = 0, 1, . . . , y n+1 = x n + dy n (4) When considering the two functional forms of h, namely exponential and rational, the above system reads: and respectively The equilibrium points of system (4) satisfy the equations andȳ =x + dȳ.
, by Then the following lemma holds.
The standard technique for analyzing the behavior of solutions near the equilibrium point as well as for exploring the dynamics under general conditions is the linearized stability analysis method, which is summarized for example in Kulenovic and Merino [14]. The analysis is based on the location of the eigenvalues of the Jacobian evaluated at the positive equilibrium (x,ȳ) relative to the unit circle. The characteristic equation reads Theorem 1. Assume that 0 < a, d < 1, α ≥ 0, and β > 0. Then, the following statements hold: of system (4) is unstable. More precisely it is a repeller.
Proof. The proof follows by checking the conditions in Theorem 2.12, p 99 in [14].
(i) Necessary and sufficient conditions for all roots of characteristic equation (10) to lie inside the unit disk are obtained if and only if Since equation (13) is never satisfied we obtain the conclusion of (ii). (iii) Verifying the conditions of Theorem 2.12, p 99 in [14] we obtain that (x,ȳ) is never a saddle point. The unstable equilibrium (x,ȳ) is a repeller if and only (15) is always true, we get the conclusion of (iii).
The following lemma gives a sufficient and easily verifiable condition for the equilibrium of system (5) to be locally asymptotically stable.
is satisfied then the positive equilibrium of system (5) is locally asymptotically stable.
Proof. From the above we have that a necessary and sufficient condition for equilibrium to be a sink is Also, remember that the equilibrium pointȳ satisfies the equation which, upon replacement in (17) produces the equivalent constraint The left-hand side of the last equation, as a function ofȳ, is seen to be monotonically decreasing. Thus, the constraint (17) is always satisfied provided that ad + aβ 1 − a + αβ < 1, and the lemma is proved.
Since we are exploring the question of what is the relative ability to effect community and individual action, we test various scenarios of condition (C1) to hold true. Because both a and d are less than 1, one possible scenario where the condition (C1) is satisfied is if α and β are very small. The parameter α is small if the impact of community consciousness on the removal of breeding sites is larger than the impact of consciousness of individual creating new breeding sites. Regarding local asymptotic stability of rational model (6), conditions in Theorem 1 translates into the following lemma.
Proof. Necessary and sufficient conditions for the positive equilibrium of system (6) to be a sink reduce to the following inequality While the left-hand side is always true, the right hand side of the above equation is equivalent to From the second equilibrium equation, we havex = (1 − d)ỹ and thus . Substituting in relation (18) we obtain Fixing y and maximizing with respect to a and d (keep in mind that the expression is linear in a and d), the maximum occurs at the extremum points. Therefore we have to consider four cases: 1) a = d = 0. The inequality (19) can be easily seen to be true by finding the common denominator. This completes the proof.
A biological interpretation of our result is worth mentioning here. As opposed to the exponential model, if the impact of consciousness decreases the formation of survival of breeding sites only gradually, then the system reaches a stable equilibrium determined by the parameters and without spontaneous outbreaks.
2. Attracting and Invariant Intervals. Global Stability Results. In this section, we obtain tight attracting and invariant rectangles for our models, intervals that are crucial in proving global attractivity results. In particular, boundedness and persistence of the systems follows. In the sequel we work with the general rescaled system (4) and recall that h is a continuously differentiable and decreasing function with h(0) = 1 and lim y→∞ h(y) = 0. Once again, the parameters satisfy 0 < a < 1, 0 < d < 1, p > 0, c > 0, α ≥ 0 and β > 0. Of importance in our presentation is the function f as defined by equation (9). The function f (y) also has the property that f (y) < 0 and, therefore, is monotonically decreasing. The and f (y) > lim y→∞ f (y) = 0 show that the codomain of f is well-defined. Then f 2k+1 (y) is decreasing, whereas f 2k (y) is increasing. 3. All fixed points of f 2 = (f • f )(y) that are different fromȳ occur in distinct pairs of the form (y 0 , f (y 0 )), with y 0 <ȳ and f (y 0 ) >ȳ. 4. The sequence {f 2k (0)} k≥0 is strictly increasing, whereas {f 2k+1 (0)} k≥0 is strictly decreasing. The two sequences converge to y m and y M = f (y m ), which are defined as the smallest and the largest fixed points of f 2 (y), respectively. Proof.

The inequality f
From the equation x n+1 = ax n h(y n ) + βh(αy n ) and the fact that h is decreasing, we have that The quantity in the curly bracket is zero by the definition of f. Indeed, which by expansion produces By a similar line of thought, we obtain which, because h is decreasing and f 2k−1 (0) > f 2k+1 (0), leads to Finally, I ∞ = ∩ k≥0 I k is invariant because it is an intersection of invariant sets. Using the first equation of the system, we can easily see that x n+1 = ax n h(y n ) + βh(αy n ) ≤ ax n + β and since 0 < a < 1 we have M 1 < β 1 − a < ∞. From the equation y n+1 = x n + dy n we obtain that ∀ > 0, there is N ≥ 0 such that y n+1 ≤ M 1 + + dy n for all n ≥ N . It follows that M 2 ≤ (M 1 + )/(1 − d) for all > 0. Since > 0 is arbitrary, we have that A similar argument with the inequalities reversed yields Next we claim that sup From the preceding equality and the definition of limsup, we have that ∀ > 0, there is N ≥ 0 such that h(y n ) < h(m 2 ) + for all n ≥ N . From the equation x n+1 = ax n h(y n ) + βh(αy n ) and the fact that h is continuous and decreasing, we Since > 0 is arbitrary, we have that A similar argument with the inequalities reversed produces Combining the equations (20) and (24), we obtain whereas from equations (21) and (25), we obtain Because f is decreasing, from equation (26), we obtain f (M 2 ) ≥ f 2 (m 2 ), which, together with equation (27) produces m 2 ≥ f 2 (m 2 ). Using that f 2 (y) is increasing, an induction argument shows that m 2 ≥ f 2k (m 2 ), whereas the fact that f 2k (y) is increasing produces m 2 ≥ f 2k (m 2 ) ≥ f 2k (0). Taking  The next corollary is a consequence of the previous lemmas.
Proof. It follows easily from Lemma 5 and Lemma 4 that lim inf n→∞ x n > 0 and lim inf n→∞ y n > 0.
According to previous lemmas, the key to obtaining global attractivity results for both systems is to establish when the solution of (f •f )(y) = y is unique. Therefore, we have the following general corollary.

Corollary 2.
If (f • f )(y) = y has a unique fixed point then the positive equilibrium is a global attractor.

Remark 1.
Even if the functions f and g are monotone in their arguments, we could not apply monotonicity theorems like the ones given in Kulenovic and Nurkanovic [15] (see also Kulenovic and Merino [14], pp 151-152 or Grove and Ladas [10]) because those type of theorems require a compact domain and codomain for our functions. It is true that we have demonstrated that the problem can be restricted to compact intervals by Lemma 4 and Lemma 5. However the proof automatically leads to the conditions appearing in the monotonicity type theorems for systems, theorems that are due to Kulenovic and Nurkanovic in [15].
A sufficient condition for the positive equilibrium of the rational system (6) to be globally asymptotically stable is given by the following theorem: Theorem 2. Assume 0 < a < 1, 0 < d < 1, α ≥ 0 and β > 0. If aα < 1 then the positive equilibrium of system (6) is a global attractor.
Proof. We are looking for the zeros of the function g(y) and h(y) = 1 1 + y . We want conditions in terms of parameters such that f (f (y)) has a unique fixed point. Then, according to Corollary 2 the result follows. Consider f (f (y)) = y, which is equivalent to Replacing f (y) in the above, we have (30) On the other hand, according to the Descartes rule of signs, the cubic polynomial has at most one real positive root. Using equation (30) we conclude that the only root of the cubic polynomial isỹ. For the quadratic term, we observe that if 1 − αa > 0, then there are no sign changes of the coefficients and, therefore, there are no real positive roots. In conclusion, if aα < 1, then f (f (y)) = y has a unique fixed point and the theorem is proved.
Local asymptotic stability of (x,ỹ) of rational system holds for a bigger parameter region than that of found analytically in Theorem 2. Combining Lemma 3 with Theorem 2 we obtain that the positive equilibrium of the rational model is globally asymptotically stable whenever aα < 1 i.e aq p < 1. In particular, it follows that we also have global asymptotic stability whenever α ≤ 1, which in terms of the original parameters is equivalent to q ≤ p, meaning that the effect of community consciousness is stronger than the individual one. In the case when aα ≥ 1, numerical simulations (see Figure 1) depict a 'spiraling in' behavior and convergence to equilibrium for system (6). In fact, simulations performed on rational system (6) indicate that the locally asymptotically stable region given in Lemma 3 is a globally asymptotically stable region. Also, computer simulations show that the behaviors described in Figures 1-3 do not appear to depend upon initial conditions. The exponential model (5) shows more complex behavior. Outside the local stabil-  ity region we encountered periodic solutions ( Figure 2) and quasi-periodic solutions ( Figure 3). It is possible that oscillations under some parameter values, with periods longer than the incubation time of the pathogen within the host, might be valuable in interrupting transmission of a mosquito borne disease. Perhaps an intervention program that focuses on this characteristic of the system might decrease the prevalence of a mosquito borne disease. Numerical simulations have also shown that system (5) supports chaotic dynamics (Figure 4).
where γ = bqc > 0. The system (32) possesses a unique positive equilibrium satisfies all the properties of Proposition 1. By applying the same techniques, results, and corollaries of section 2, the global attractivity region of the equilibrium point of system (32) is given by the condition that the equation (f •f )(y) = y has a unique fixed point. We look at the two special forms of function h, namely h(y) = h 1 (y) = 1 1 + y (rational) and h(y) = h 2 (y) = e −y (exponential).
For the rational case we have the following theorem: Theorem 3. Assume 0 < a < 1, 0 < d < 1, and γ > 0. Then the positive equilibrium of system (32) is a global attractor.
Proof. We claim that the equation = y has a unique solution for the values of the parameters, as in the hypothesis. Replacingf (y) in the above and simplifying we obtain that For the exponential case, the following theorem holds: , then the positive equilibrium of system (32) is a global attractor.
Proof. The equation (f •f )(y) = y is equivalent to . Consider the functiong(y) = Ae − A e y − y. We have that g(0) > 0, and lim y→∞ g(y) = −∞ and therefore there exists at least one positive root ofg(y) = 0. On the other hand, because g (y) = A 2 e −Ae −y −y − 1 < 0 whenever A < 1, we obtain that g is decreasing and thereforeg(y) = 0 has exactly one root on the interval (0, ∞).
Remark 2. : The condition γ < (1−a)(1−d) translates into bqc < (1−a)(1−d) and the previous theorem states that global attractivity is obtained when the product bqc is very small. One possible scenario when this can happen is when b and c are big and q < 1 bc , i.e., q is very small, meaning that individual consciousness has a small impact on the number of breeding sites.
4. Concluding Remarks. The described models are an attempt to represent community involvement in mosquito control. Mosquito control through habitat management, for the purpose of protecting the population from Dengue fever, has been a common practice in tropical and subtropical countries. This is implemented by educational campaigns raising individual and community consciousness, as captured by the two different nonlinear systems of difference equations presented in this work. In the two models community and individual awareness are in a negative feedback with mosquito breeding sites, allowing for rich dynamics. In both the exponential and rational models, there is a parameter region in which the system is globally asymptotically stable, but this region is larger in the rational model, where the impact of the awareness on the system behavior is weaker.
From the point of view of the administrative authorities, the dynamics must be brought in a region in which the results are predictable (that is, globally stable) and for which the number of breeding sites is consistently low. The first requirement is motivated not only by the financial constraints but also by the planning strategies, which normally desire definite and measurable outcomes. The tendency of the exponential system to exhibit chaotic behavior shows that a high impact of the community and individual awareness on the dynamics of the system is not actually desirable, since the outcomes may not depend continuously on the incremental changes the administrator is likely to implement. This lack of consistency in actual results is likely to stress the budget and wear down the very social awareness one is striving to build up.
The global stability results we have been able to demonstrate suggest that a strategy with predictable outcomes is obtained by less aggressive interventions on the public awareness, as modeled by the rational system. Sure, "less aggressive" does not mean that the number of breeding sites will become smaller, but asymptotically stable, only. The actual number of breeding sites can be diminished by controlling the parameters defining the unique global equilibrium. The less aggressive policy has to be maintained for extended periods of time, in order for the equilibrium to be attained. Also, a public campaign must always be matched as magnitude by actual efforts of the administrator to remove a portion of the breeding sites. This is so because our global stability results for α < 1 show that the effect of community consciousness is stronger than the effect of individual consciousness in ensuring predictable outcomes. On the other hand, the number of breeding sites does decrease as α increases closer to 1. This is apparent in Figure 4 even for the exponential model (if restricted to its small region of global stability). As such, public education should not be given up altogether, especially if it is a cheaper alternative compared to the actual removal costs.