eu PREDATOR-PREY MODEL FOR ANALYSIS OF AEDES AEGYPTI POPULATION DYNAMICS IN CALI , COLOMBIA

In this paper, we introduce a predator-prey model to analyze the population dynamics of the Aedes aegypti mosquito which is fairly blamed for the transmission of dengue in Cali, Colombia. The model describes an interaction of an aquatic predacious species deliberately introduced into local mosquito breeding habitat, where the Aedes aegypti immature stages (eggs, larvae, and pupae) are already present and constitute the prey population. The model also accounts for local population of adult female mosquitoes (or mature stage) emerging from the breeding site. This population is considered as a target for reduction by deploying an adequate predacious species since only female mosquitoes are held responsible for transmission of dengue and other vector-borne diseases. Having analyzed the model, we have derived the mosquito survival threshold with predation as a function of predator’s biological characteristics. The model’s parameters were adjusted to the average seasonal temperatures of Cali, Colombia and explicit conditions for biological characteristics of prospective efficient predators were established. Numerical simulation with introduction of an efficient predacious species in local mosquito breeding habitats revealed the possibility of eventual mosquito extinction in such localities. Finally, some particular biological species were proposed as potential candidates for efficient predators. Received: April 29, 2015 c © 2015 Academic Publications, Ltd. url: www.acadpubl.eu Correspondence author 562 J.H. Arias, H.J. Martinez, L.S. Sepulveda, O. Vasilieva AMS Subject Classification: 92D25, 93D20, 34A34


Introduction
Aedes aegypti is an invasive mosquito species (family Culicidae, order Diptera) that has colonized all tropical and sub-tropical regions worldwide where this species is fairly blamed for transmission of vector-borne diseases [4].After mating, the female mosquitoes need human blood meals in order to mature their eggs, and can easily transmit the pathogen or virus between people they feed on.
Aedes aegypti is strongly anthropophilic, peridomestic and opportunistic mosquito species; its breeding sites and lifetime habitats can be found in and around habitable houses, workplaces and public areas with constant people flow or congregation.Therefore, the presence and abundance of Aedes aegypti populations in many tropical countries (and notably in Colombia) are strongly correlated with dengue infections, and many Colombian cities are considered hyperendemic with regards to dengue morbidity [23,29].One of them is the city of Santiago de Cali (referred to as Cali from here on), as it is pointed out in [14].Figure 1.1 shows the dengue incidences (i.e.registered cases) collected by the Municipal Public Healths Secretariat of Cali during 2001-2013, and discloses an alarming trend in propagation of this disease.
In absence of approved vaccine (see, e.g.[15]), the most effective intervention to reduce dengue infections consist in targeting the mosquito population through the use of insecticides and larvicides.This method has been routinely deployed all over Colombia for more than five decades [38].In Cali, the principal breeding sites of Aedes aegypti are the rainwater catch basins located along street in residential and commercial areas, and they are constantly treated with chemical substances [27] 1 .Nonetheless, the dengue morbidity continues to increase.On the one hand, this unfortunate consequence can be attributed to insecticide/larvicide resistance acquired by mosquitoes when certain chemical substances are being applied during considerable time [11,28,39].On the other hand, the abundance of Aedes aegypti populations that naturally causes higher dengue incidence is typically correlated with ongoing climatic changes (e.g., higher temperatures, humidity, and rainfalls) that provide more favorable condition for mosquito breeding and reproduction [1,16,41,46].
The actual situation in Cali (see again Figure 1.1), where chemical eradication actions have never been suspended during last decade, clearly illustrates that a sole use of insecticides and larvicides is not sufficient to bring dengue infections under control.Global analysis of overall dengue propagation also corroborates the tendency observed in Cali [12,13].
Despite of the resistance threats, environmental concerns, and possible toxicity to non-target organisms (including human population), chemical control still remains the pillar of vector control programs in many countries (see, e.g.[5,10,22] and references therein).However, the World Health Organization (WHO) and the Pan-American Health Organization (PAHO) jointly promote and encourage its combination with other non-chemical vector control strategies2 .Among them, we would like to emphasize the biological control strategy, i. e. the use of living organisms that interact with mosquito populations by means of predation, competition or parasitism, as well as genetic modifications of mosquitoes seeking to reduce their longevity and/or fecundity or, otherwise, to disable them for disease transmission [5,22,34].
Biological control is generally regarded as an environmentally friendly method to mitigate risks and negative impacts in ecology, agriculture, or public health caused by global factors.In this sense, it is worthwhile to study the Aedes aegypti population dynamics in the presence of biological predators of the imma-ture stages since they are of low cost and do not contaminate the environment.The allocation of proper predacious species in targetted Aedes aegypti breeding sites can be performed by a single release (also regarded as inoculative) or by multiple inundative releases as a complementary measure to traditional chemical control over the mosquito population.
This paper proposes a mathematical model of predator-prey type that describes the population dynamics of two mosquitos stages -immature or aquatic (comprising the eggs, larvae, and pupae) and mature or aerial (encompassing only female adult mosquitoes) 3 .
This model has local nature meaning that it focuses on the population dynamics inside and around a generic breeding place (referred to as a "basin" from here on) located in some public area that accounts for constant people flow or congregation.Examples of such "basins" are widely spread across many tropical cities, including Cali, and form an essential part of their conceptual architectonic design.In particular, Cali has various shopping and recreational centers, parklands and clubs, residential complexes, open-air restaurants, school and university campuses decorated with ornamental concrete and tiled bowls, basins, pools and other cavities filled with stagnant water and surrounded by rich vegetation.From the aesthetical standpoint, these landscapes are spectacular; however, they constitute perfect sites for mosquito breeding and thus indirectly contribute to transmission of vector-borne diseases.It should be pointed out that property owners and/or managing authorities (at least, in Cali) are strictly opposed to the use of chemical substances (insecticides and larvicides) within these premises since such measures would imply a temporal closure of their establishments, resulting in a negative impact on its images and scaring off some potential visitors or customers.However, local authorities may easier agree to environmentally acceptable method of biological control (use of natural predators within the basins located in their territories) and thus make their premises safer for all visitors.
The model (presented in Section 2) is purposely designed in order to provide valid arguments in favor of biological control by using natural predators of the Aedes aegypti immature stages in multiple locations around Cali where no other control measures are being applied.Thorough analysis of the model (performed in Sections 3 and 4) allowed to derive essential conditions to be imposed on biological characteristics of an appropriate predacious species capable to reduce the local Aedes aegypti populations.Our rationale is straightforward: smaller populations of Aedes aegypti effectively present in the crowded public areas should be reflected in lower dengue incidence since this mosquito species is fairly held responsible for the disease transmission.Section 5 provides estimations of the model parameters with regards to three average temperature of Cali (corresponding to the cool, warm, and hot seasons) and indicates, in numerical terms, what biological characteristics of predacious species are required for efficient biological control of Aedes aegypti in Cali.
Numerical simulations (Section 6) performed for different sets of model parameters help to envisage the benefits of using an appropriate predator (whose biological features comply with imposed conditions) in generic basins located across Cali.Further discussion provides some ideas regarding a particular predacious species to be employed as a biological control agent in Cali.Finally, Section 7 indicates how can the outcomes of the paper be used in the practical context in order to reduce the Aedes aegypti population actually present in Cali, Colombia.

Description of the Model
Before formulating the model, one should clearly understand the reproduction mechanism and life cycle of the mosquito populations.Aedes aegypti females lay their eggs on sidewalls (just above water surface) of natural stagnant waterbodies (ponds, swamps, tree holes, leaf axils, fallen fruit husks, etc.) and water-storage containers used in and around households (flower vases, discarded tires, buckets, cisterns, ornamental basins, etc.).Eggs may remain viable up to several months waiting for the water to cover them and reach an adequate temperature for their further hatching into larvae (in tropical climate it usually occurs within 2-7 days).The larvae (comprising 4 instars) are transformed into pupae within 4-6 days and the pupae eventually become adult mosquitoes in 1-3 days.On the second day after emergence, the adult mosquitoes are able to copulate and the eggs are matured for 2-4 days.Each female mosquito may have up to 6 gonotrophic cycles throughout her lifespan, and multiple ovipositions may occur during the same gonotrophic cycle.Figure 2.1 shows the Aedes aegypti life cycle 4 and further details regarding entomology, bionomics, and life cycle of Aedes aegypti are meticulously described in [7,32,44].
Following the framework proposed in [46] for population growth of the Aedes aegypti, we consider two mosquito stages -immature or aquatic (compris- ing eggs, larvae, and pupae) and mature or aerial (consisting of adult female mosquitoes).Here we suppose that there are enough male mosquitoes in the surroundings of a generic basin (breeding site) to mate with females.Therefore, the male subpopulation is ignored in the modeling, assuming that all female mosquitoes are able to copulated prolifically.Additionally, we consider a predator-prey interaction between the immature stages (prey) and some aquatic predacious species deliberately placed in the basin (breeding site).Thus, we have three state variables: m q (t) -number of immature stages (eggs, larvae, pupae) at the moment t, m a (t) -number of adult female mosquitoes at the moment t, x(t) -number of natural predators of immature stages at the moment t.In mathematical terms, the resulting model is stated as follows: ) ) with initial conditions m q (0) = m q 0 , m a (0) = m a0 , x(0) = x 0 , whose positive (constant) parameters are described in the Table 2.1.Eq. (2.1a) basically states that immature stages m q (t) (eggs, larvae, and pupae) proliferate due to oviposition rate of the adult female mosquitoes (ε), the fraction of the eggs hatched (k), the number of female adult mosquitoes (m a ) capable to lay eggs, and the remaining carrying capacity (1 − m q /K q ).The volume of immature stages decays proportionally to their transition into mature stage (νm q ), due to natural mortality (δ q m q ), and as a result of predation (ϕxm q ) which is expressed using the "law of mass action".We presume that predation of immature mosquito stages is carried out regardless of the prey gender; in other words, both male and female immature stages are being predated.Therefore, a fraction f of adult female mosquitoes is considered for m a (t) dynamics rather than for m q (t), in contrast to original model developed in [46].Similarly, Eq. (2.1b) indicates that the number of adult female mosquitoes m a (t) increases with transition from immature to mature stages (νm q ) where only a fraction of female mosquitoes (f ) is taken into account.The population of m a decreases due to its natural mortality (δ a m a ).
A single specie of natural predator x is placed in the basin in order to reduce the population of adult female mosquitoes m a which are potentially capable to transmit the dengue virus.Thus, the ultimate role of x consists in reducing the propagation of this infectious disease.Eq. (2.1c) describes logistic growth of such predator and states that the predator's diet is not limited to the mosquito's immature stages.The latter implies that x has alternative food sources such as algae, ciliates, other insects and their immature stages that may be present at the same habitat.Therefore, we consider an additional increase ψm q in predator's intrinsic growth attributed to its interaction with this particular prey.
Remark 1.The authors of [19] had made an attempt to include an acting predator into 2-stage mosquito life cycle model initially proposed in [46].However, in their version the predator-prey interaction was assumed unilateral, meaning that the predator population was supposed to grow independently of the prey population.Under this assumption, the differential equation for predator population growth (identical to (2.1c) for ψ = 0) was simply decoupled from the other two, and a two-state system was studied.
It is also supposed that the predacious species occupies the total volume of the water in the basin, while the habitat of m q is limited to the basin's walls (for eggs) and to the water surface (for larvae and pupae).Therefore, the carrying capacity of the predator, K x , is usually greater than the carrying capacity of the immature stages, K q , if the individuals of x and m q have comparable sizes.
Additionally, model (2.1) is of local nature in the sense that it describes the growth of Aedes aegypti population around some generic basin, since the interaction between m q and x takes place within such basin and adult female mosquitoes usually do not fly far from the basin (or container) where they developed as larvae and pupae 5 .
It is worthwhile to note that the region of biological interest, that is, is positively invariant.The latter means that any trajectory of the system (2.1) that starts in Ω will remain in Ω for all t ≥ 0. A formal proof of this statement can be found in the Appendix A.

Stability Analysis
To analyze stability of the dynamic system (2.1), we calculate its points of equilibrium and then formulate a criterion that permits to establish local stability properties of all equilibria.Clearly, the points of equilibrium are solutions of the following algebraic system of nonlinear equations: From Eq. (3.1b), we have that By replacing (3.2) in the Eq.(3.1a), we can reduce the system (3.1) to Eq. (3.3b) gives two options, namely: and leads us to two possible cases: x = 0 (predator's extinction or absence) and x = K x (predator's persistence).These two cases are briefly revised below, and their analysis results in two useful thresholds for mosquito's survival that rule upon stability of the underlying equilibria.
• Option 1: x = 0 at the point of equilibrium.In this case, we have from Eq. (3.3a) that The last expression will bear biological sense (that is, m q > 0), only if where Q can be viewed as mosquito's survival threshold without predation6 that represents an average number of female mosquitoes produced by one female mosquito in the absence of predator of the immature stages (eggs, larvae, pupae).
Using this notation, the component m q can be expressed as Thus, in virtue of (3.4) and (3.6), we have two options: - Combining them with the absence of predator (x = 0), we obtain two possible equilibria • Option 2: x = K x at the point of equilibrium.In this case, Eq. (3.3a) implies that The last expression will bear biological sense (that is, where Q 0 can be viewed as mosquito's survival threshold with predation that indicates an average number of female mosquitoes produced by one female mosquito in the presence of predator of the immature stages.One may also observe that (εkf ) δ a in (3.9) is an average number of larvae that turn into female mosquitoes, and ν (ν + δ q + ϕK x ) is the probability that such larvae survive, despite of predation, and become adults.
Remark 2. It is clear that since both ϕ and K x are positive.
Using notation (3.9), the component m q can be expressed as Thus, in virtue of (3.8) and (3.11), we have two options: - Combining them with the presence of predator (x = K x ), we obtain two additional equilibria Observe that, among all four points of equilibrium, only E 4 admits coexistence between the prey m q and predator x.
Existence and stability properties of the equilibria (3.7), (3.12) naturally depend on the thresholds values Q, Q 0 and these properties can be summarized in the following theorem.
Theorem 1.In view of the inequality (3.10) we have three possible options for the system (2.1), namely: and E 3 is an attractor (locally asymptotically stable equilibrium).Note that E 2 and E 4 have no biological sense in this case.

If
, and E 3 is an attractor (locally asymptotically stable equilibrium).Note that E 4 has no biological sense in this case.

If
), E 2 and E 3 are saddle points (conditionally stable equilibria), and E 4 is an attractor (locally asymptotically stable equilibrium).
Formal proof of the above theorem is provided in Appendix B. Theorem 1 clearly states that only two equilibria (E 3 and E 4 ) are potential attractors of the predator-prey system (2.1) and their local stability directly depends on the value of Q 0 .Thus, under initial presence of predator (that is, x(0) > 0) together with at least one stage of mosquito (either m q (0) > 0 or m a (0) > 0), the trajectories of the system (2.1) will converge either to E 3 (mosquito-free equilibrium) or towards E 4 (coexistence).
Apparently, the value of Q 0 (cf.(3.10)) depend on biological characteristics of the predator -its carrying capacity K x and efficiency of predation ϕ.
If one female mosquito produces, in average, less than one female mosquito due to timely predation of the immature stages (that is, Q 0 < 1), then the underlying predator should be referred to as efficient and its presence in the basin will eventually lead to local extinction of the Aedes aegypti population around the basin.
Otherwise, if one female mosquito produces, in average, more than one female mosquito regardless of predation of the immature stages (that is, Q 0 > 1) the predator is regarded as inefficient and its presence in the basin will eventually lead to coexistence between prey and predator thus guaranteeing local persistence of the Aedes aegypti population.
Therefore, by choosing a predator with appropriate biological characteristics (K x and ϕ) one may considerably reduce the value of Q 0 and favorably change the asymptotic behavior of the system (2.1).

Sensitivity Analysis
A formal sensitivity analysis will help us to learn the impact of each parameter upon the global value of the mosquito's survival threshold in the presence of predator .
All eight parameters of Q 0 can be eventually varied by external forces.For example, according to [46] the oviposition (ε), egg-hatching (k), transition (ν), and natural mortality (δ a , δ q ) rates are sensitive to climate changes such as temperature fluctuations, droughts, scarce or excessive precipitations, etc. Mortality rates of immature stages (δ q ) and adults mosquitoes (δ a ) can be also controlled using larvicides and insecticides, respectively 7 .By choosing a predacious species with adequate biological characteristics ϕ and K x , one may achieve a noticeable reduction in the value of Q 0 .
Mathematically, Q 0 can be viewed as a function of eight variables, and its relative change with respect to each particular variable can be expressed by the According to the above expressions, we have (ceteris paribus) that Q 0 is an increasing function of ε, k, f, and ν, while Q 0 is a decreasing function of δ a , δ q , ϕ, and K x .It is also useful to employ the concept of normalized forward sensitivity index (NFSI).According to [6], the normalized forward sensitivity index of a variable, V , that depends differentially on a parameter, p, is defined as This explicit formula permits to compare the incidence of each parameter over the global value of mosquito's survival threshold in the presence of predator Q 0 .Table 4.1 presents the NFSI values with respect to each parameter.It shows that, for example, if the oviposition rate, ε, is increased (or decreased)  by 10%, then, Q 0 increases (or decreases) also by 10%.We can observe that any possible changes in both predator's parameters K x , ϕ will have no impact on oviposition rate ε, fraction of hatched eggs k, fraction of female mosquitoes f , and natural mortality rate of the adult mosquitoes δ a since their corresponding indexes remain constant.
On the other hand, the normalized forward sensitivity indexes (NFSIs) Υ ν , Υ δq , and Υ Kx = Υ ϕ explicitly depend on the product ϕK x (see Table 4.1), and Figure 4.1 provides their absolute-value plots as real functions of the product ϕK x (note that Υ ϕ = Υ Kx are displayed by the same curve).In addition, the change in these four NFSIs is stronger in the lower range (that is, between 0 and 0.5), while in the upper range of the product ϕK x (between 0.5 and 1), all curves have a more pronounced asymptotic behavior.
It is also convenient to re-write the mosquito's survival threshold with predation, Q 0 , in the following way: , where A = εkf ν δ a and B = ν +δ q (4.1) are constants defined by the mosquito's entomological characteristics.In line with [46], these characteristics (as well as Q) are temperature-dependent and higher temperatures provide better conditions for mosquito's survival and re- production.As a result, the value of Q 0 should also depend on ambient temperature 8 .
Therefore, seasonal variations of the entomological parameters conforming A and B in Eq. (4.1) should be taken into account when choosing a predator with appropriate biological characteristics ϕ and K x .

Aedes Aegypti Population Dynamics in Cali, Colombia
Cali belongs to hyperendemic regions of Colombia with regards to dengue morbidity and the Aedes aegypti population is widely spread around its urban and semi-urban areas [14,23,29].Cali is located 1 000 meters over sea level in the valley (called Valle del Cauca) close to the western Andean mountain chain and occupies about 560 square kilometers.According to 2005-2020/DANE population projection (see more details in [8]), Cali has about 2.4 million inhabitants.
Due to its proximity to the equator, there are no major seasonal variations.Over the course of a year (see Figure 5.1), the average temperature may vary from 19 o C (cool periods with abundant rains, 1 to 3 months per year) to 29.8 o C (hot periods with possible drought, up to 4 months per year) but it typically remains at about 23.9 o C during most part of the year.
Cali has flourishing vegetation all year around and its original concept of urban design stems from the landscapes where ornamental concrete or tiled water basins are surrounded by rich vegetation.From recreational and aesthetical points of view, these landscapes are spectacular.However, artificial basin filled with stagnant water are perfect sites for mosquito breeding, and such basins are widely spread across the city.A great part of them are located in crowded public areas (such as shopping and recreational centers, residential complexes, parklands, schoolyards, university campuses, etc.) where all visitors are in risk to be bitten by mosquitoes and possibly become infected with dengue.
Our model (2.1) describes local population dynamics of the Aedes aegypti mosquito in and around one single basin.Therefore, this model can actually help us to envisage the population dynamics of Aedes aegypti in and around a "generic basin" located in Cali, Colombia.For that purpose, the model's entomological parameters should be properly adjusted to seasonal average temperatures of Cali.
The authors of [46] held temperature-controlled experiments and thoroughly described the procedure for estimation of temperature-dependent entomological parameters of the Aedes aegypti mosquitoes which is based on polynomial fitting by the linear least squares method.This technique was applied to three average temperatures of Cali under two reasonable assumptions: (a) almost all the eggs of Aedes aegypti hatch into larvae (that is, k = 0.9); (b) half of them are developed as females (f = 0.5).
The outcome data are presented in Table 5.1 and underlying technicalities (omitted here) can be consulted in [2].Table 5.1 also displays the mosquito's survival threshold without predation Q that naturally increases for higher temperatures and always remains well above 1 for all considered temperatures.This implies that, under no vector control actions, the Aedes aegypti population will never become locally extinct in Cali, Colombia.This situation can be remedied by a deliberate and simultaneous introduction of some efficient predacious species in ornamental basins (that are vastly spread around Cali) in order to decelerate the process of mosquito proliferation.
In mathematical terminology, one should find a predacious species, x(t), whose biological characteristics (K x and ϕ) would comply with condition Q 0 < 1 preferably for all average temperatures considered.Here Q 0 can be viewed as a function of two variables ϕ and K x , while other parameters may remain constant within each season (see Eq. (4.1)).
Figure 5.2(a) displays the level set of Q 0 (ϕ, K x ) for warm-season set of parameters (see Table 5.1, T =23.9 o C) 9 .This graph has an obvious interpretation: larger carrying capacity of the predator K x allows for lower effectiveness of predation ϕ and, similarly, smaller carrying capacity K x requires higher effectiveness of predation ϕ in order to maintain constant the value of Q 0 .Remark 3. In economics terminology, these two parameters K x and ϕ are interchangeable (or, more precisely, intersubstitutable) since for Q 0 in the form (4.1), the marginal rate of substitution (MRS) of ϕ for K x is calculated as and indicates that by giving up (K x /ϕ) units of the carrying capacity of predator K x , one can get one more unit of predation effectiveness ϕ while holding constant the value of Q 0 .
Applying the MRS concept, we can observe that during warm season (see Figure 5.2(a)) a 10% drop in predation rate ϕ from (0.3 to 0.2) will require to double the carrying capacity (from 30 to 60 units) in order to hold the value of Q 0 slightly below 1.
Since ϕ and K x are intersubstitutable, their product, ϕK x , can be considered as a single variable.Figure 5.2(b) helps to visualize Q 0 as a real function of the product ϕK x with regards to cool, warm, and hot average seasonal temperatures in Cali, Colombia.In this graph, we can see what values of the product ϕK x are capable to maintain Q 0 below 1 and observe that these values should significantly increase for higher temperature.Indeed, we can see that for coolseason average temperature (19 o C), the product ϕK x should be greater than 1.6; for warm-season average temperature (23.9 o C), ϕK x should be greater than 5.45, and for hot-season average temperature (29.8 o C), it must exceed 18.31.
Following this idea, one may define a "target value" of the product ϕK x and select an "appropriate predator", i.e., a particular predacious species whose biological characteristics comply with such target value.For example, setting ϕK x > 15 as target, and after infesting our "generic basin" with a predacious species that satisfies this condition, we may expect that local Aedes aegypti population (both immature and mature stages in and around the basin) will steadily decline during cool and warm seasons.However, during hot season a coexistence between prey and predator can also be reached.

Numerical Simulations and Discussion
In this section, we present the numerical solutions of the system (2.1) which were obtained using the MATLAB ode45 routine.Numerical values of mosquito's entomological parameters are taken from the Table 5.1 for different average temperatures and assuming that k = 0.9, f = 0.5 (i.e., 90% of all eggs hatch into larvae and the half of them are developed then as female mosquitoes).The carrying capacities of immature stage and predator have been normalized to K q = 1 and K x = 10, respectively 10 , and initial conditions for mosquito have been assigned as fractions, that is, m q (0) = m q 0 = 0.2, m a (0) = m a0 = 0.3.(6.1) By setting K x = 10, another predator's parameter ϕ (predation rate characterizing the efficiency of predation) will define two alternative types of possible predators.Namely, (1) if for chosen value of ϕ the condition Q 0 > 1 holds, then the corresponding predator enables predator-prey coexistence and, therefore, it should be classified as incapable or inert; (2) if for chosen value of ϕ the condition Q 0 < 1 holds, then the corresponding predator is capable to induce mosquito's local extinction and, therefore, it should be classified as efficient or active.
In theory, local extinction of Aedes aegypti can be achieved by introducing an efficient predator whose biological characteristics ϕ, K x comply with the following condition: It was supposed that a generic predator population grows a bit faster due to timely predation of the Aedes aegypti immature stages (see Eq. (2.1c)) and its additional growth rate is expressed by ψm q .It is worthwhile to note that ψ > 0 has no incidence on the mosquito's survival threshold with predation, Q 0 , so this parameter can be fixed for further simulations at some reasonable level, e.g., ψ = 0.1.Roughly speaking, any increase in ψ will slightly speed up the growth of predator population towards its carrying capacity (due to extra food available) without altering the stability features of prey population.
Finding a predator whose biological characteristics would comply with (6.2) for all three average temperatures seems as a challenging task.In that sense, we consider more plausible to find a particular predacious species whose biological features will fulfill the condition (6.2) for warm-season temperature (and, naturally, for cool season too) while this condition will be violated during hot season.
Predator's population dynamics is quite predictable (logistic growth for any x(0) = x 0 > 0) and is of no particular interest here; therefore, in all consequent charts the trajectories of x(t) are suppressed.Our simulations are principally focused on comparison of two possible scenarios: local population dynamics of immature stages m q (t) and adult female mosquitoes m a (t) in and around some "generic basin" (located in Cali, Colombia) with and without predation of immature stages.
Supposing little seasonal temperature variations (i.e. using the parameter set of Table 5.1 for T =23.9 o C) together with initial conditions (6.1), we generate two sets of trajectories m q (t) and m a (t) of the system (2.1) for x(0) = x 0 = 0 (absence of predator) and for x(0) = x 0 = 3 > 0 (presence of predator).First set of trajectories is obtained for an inert predator whose biological features do not comply with condition (6.2) at all (that is, for all seasonal temperature variations).The second set of trajectories is obtained for relatively active predator whose biological features fulfill the condition (6.2) during cool and warm seasons.Figure 6.1 displays the results where the trajectories for m q (t) and m a (t) (green and red curves, respectively) are drawn by solid lines in the absence of predator and by dashed lines in the presence of predator.The contrast between two charts is striking.Figure 6.1(a) shows that very moderate reduction of both Aedes aegypti populations can be gained when an inert predator is introduced into basin.Alternatively, Figure 6.1(b) reveals quite notable reduction of both m q (t) and m a (t) as time goes by.
It is easy to foresee that the same predator should perform even better during cooler periods, that is, when air temperature drops from T = 23.give a reasonable answer to this question we have performed another series of simulations and their results appear in Figures 6.2 and 6.3 under the following legend: • blue-colored curves refer to cool season with average temperature T min = 19 o C; • green-colored curves refer to warm season with average temperature T = 23.9o C; • red-colored curves refer to hot season with average temperature T max = 29.8o C; • trajectories drawn by solid lines are generated in absence of predator (x(0) = x 0 = 0); • trajectories drawn by dashed lines are generated in presence of predator (x(0) = x 0 > 0).
As expected, the presence of an inert predator shows insignificant impact on both Aedes aegypti populations, m q (t) and m a (t), during warm and hot periods (see Figures 6.2(a) and 6.3(a)).However (and quite surprisingly!), the same charts suggest that even an inert predator can induce a noticeable decrease in both mosquito populations if cool temperature is sustained for some time.What possible explanation may have this phenomenon?We have not made any assumption for predator to gain additional efficiency during cool periods.Table 5.1 shows that entomological parameters are sensitive to temperature alterations.Namely, oviposition and transition rates (ε and ν) decrease with temperature drops, while mortality rates of immature stages and adult mosquitos (δ a and δ q ) increase.As a result, the mosquito survival threshold Q tends to decline for lower temperatures.However, this is also valid for trajectories m q (t) and m a (t) generated for the same set of parameters (ε, ν, δ a , δ q ) in the absence of (inert) predator (cf.blue solid-line curves in Figures 6.2(a) and 6.3(a)).If the presence of an inert predator makes almost no difference, why these trajectories, especially in Figure 6.2(a), are so different?
The correct answer stems from the Theorem 1, item 3. Obviously, Figures 6.2(a) and 6.3(a) refer to the situation when 1 < Q 0 < Q.Therefore, if x(0) = x 0 = 0 (absence of predator) then the trajectories of the system (2.1) should reach the saddle point E 2 given by Eq. (3.7b).On the other hand, if x(0) = x 0 > 0 (presence of predator) then the trajectories of the system (2.1) should reach the point of coexistence E 4 given by Eq. (3.12b).Comparing first two coordinates of E 2 and E 4 we observe that The latter explains the difference between blue curves in Figures 6.2(a) and 6.3(a).On the other hand, the presence of more efficient predator (which is referred to as "active" during cool and warm seasons) has significant impact on both Aedes aegypi populations (see Figures 6.2(b) and 6.3(b)).Figure 6.2(b) shows steep drop of Aedes aegypi immature stages m q (t) for all temperature ranges during first two week after infesting the basin with predacious species, which is followed by steady decline during cool and warm seasons.However, during hot season the situation is different.Effectively, higher temperatures facilitate smooth recovery of m q (t) thus inducing coexistence between the prey and predator.
As an outcome (see Figure 6.3(b)), the population of adult female mosquitoes m a (t) will experience steady decay during cool and warm seasons that may eventually lead to their local extinction (if appropriate temperature regimes are sustained for sufficient time).Nonetheless, the population of adult female mosquitoes will continue to grow during hot season.However, it will increase twice more slowly than in the absence of predator.Therefore, the presence of "relatively active" predator (which may become inefficient during hot season) should also have beneficial effect on overall reduction of m a (t).Summarizing all arguments presented above we arrive to the following conclusion.Since warm season generally prevails during the year, it is reasonable to expect that a "relatively active" predator (i.e., that fulfills the condition (6.2) for T = 23.9o C in context of Cali, Colombia) be capable to perform adequate control over the population of adult female mosquitoes during entire year.
Even though our arguments look solid from the mathematical point of view, there is still one pending question.Namely, what particular biological species can be considered as candidates for "relatively active" predators?
Frankly speaking, all our simulations has been performed using artificial values for ϕ and K x ; therefore, we cannot provide an unequivocal answer to this question straightforwardly.A thorough study of biological features of possible predators is required in order to decide which of them would better combat the mosquito population.Such study naturally lays beyond the scope of this paper and can be proposed for further research to professional biologists.
On the other hand, many scholars had claimed that large species of cyclopoid copepods (such as Macrocyclops, Megacyclops, and Mesacyclops) are more effective for practical mosquito control than any other invertebral predator of mosquito larvae (see [21] and references therein).Curiously, these tiny crustaceans (about 1-1,5 mm long) possess some remarkable biological characteristics, which are clearly explained in the textbook [20].In particular, it is worthwhile to mention some of them: • Cyclopoid copepods are comparable with Aedes aegypti larvae in body size and length and their natural habitat comprises the water column, water surface, and bottom sediments.Therefore, the carrying capacity of this predacious species, K x , should notably exceed the carrying capacity of the mosquito's immature stages, K q , with regards to the same generic basin.
• Cyclopoid copepods are aggressive predators and they can readily consume preys up to about twice their size.If the mosquito larvae are abundant, cyclopoid copepods eat only a small part of each larva, giving each copepod the capacity to kill 30 to 40 larvae per day, far more than they actually can eat.Therefore, their efficiency of predation, ϕ, is significantly high.
• If present in a high percentage of breeding places, cyclopoid copepods are capable to cause local eradication of mosquito by killing 95-100% of larvae.
In summary, the biological features of cyclopoid copepods fully comply with assumptions made for our model.Therefore, this particular predacious species makes a good candidate for potential predator considered in our model.

Further Considerations and Conclusions
Traditional approaches to control of vector-transmitted infections are generally focused on reducing the mosquito abundance.These methods fully rely on the use of chemical substances (insecticides and larvicides) that are spread around mosquito breading sites and habitats.Aedes aegypti mosquitoes prefer to breed in areas with stagnant water such as flower vases, uncovered barrels, discarded tires, water-storage tanks, wells, pools, artificial basins, and other water containers usually located around residential, commercial, and recreational zones.The excessive spread of chemical substances in public and domestic areas contributes to air pollution, contaminates the environment, and may provoke intoxication, allergic reactions, and respiratory disorders among human population.
In recent years, many researchers have been working on alternative vector control strategies for Aedes aegypti mosquitoes such as: 1. Genetic modification also known as the Sterile Insect Technique or "SIT" (see, e.g., [9]).
Obviously, items 1, 2, and 3 should have higher implementation costs than item 4 since they would necessarily require professional training of the personnel in charge.In practical context, last option (item 4) looks more accessible and can be carried out under relatively modest funding 12 .Additionally, numerous field and laboratory experiments have provided weighty arguments for justifying the use of cyclopoid copepods as biological agents capable to eradicate local Aedes aegypti (and other mosquito) populations.These experiments have been conducted in different natural and artificial mosquito breeding sites located in diverse regions and climatic conditions.
For example, several laboratory and field experiments had been held with cyclopoid copepods using discarded tires in subtropical climate (east-central Florida, USA) and the reported reduction in larval survival was close to 90% [31].According to [36], the experiments held in Costa Rica using bromeliad leaf axils, tree holes, and discarded tires as mosquito breeding habitats had shown the larval reduction of 79%, 90%, and 99% in tropical dry, moderate, and humid climate zones, respectively.Cyclopoid copepods were also recommended as cost-effective, environmentally acceptable and persistent agents for the sustainable control of Aedes aegypti, especially in inaccessible breeding sites such as Australian gold mines and wells [35].As already mentioned, a routine operational use of cyclopoid copepods is practised in hundreds of villages in Vietnam.The copepods are settled in all principal breeding sites of Aedes aegypti (such as communal water-storage tanks, wells, and household water containers) and the participating communities do reportedly achieve, in average, a 98% reduction of overall Aedes aegypti population (see more details in consecutive monitoring reports [17,18,24,26]).
Yet, there is evidence disclosing that populations of cyclopoid copepods may have vulnerable persistency in diverse environments.During 1999-2000, the colonies of cyclopoid copepods were introduced in some selected street catch basins of Cali, Colombia [40].The copepods had survived only in less than half of the catch basins and they maintained there low densities of Aedes aegypty larvae during the 8-month study (reported efficacy was 90.4% in larval mortality).Nonetheless, the general outcome of this experiment was not satisfactory, since the majority of copepods did not tolerate water contaminated with detergents that came from drainage of the neighborhoods, and its toxic action had killed the copepods without killing the mosquito larvae.This experiment had not fulfilled the expectation of local public health and sanitary control authorities and discouraged further use of cyclopoid copepods as potential biological agents of mosquito control in Colombia where mosquito populations are principally controlled using traditional methods (insecticides and larvicides).
On the other hand, there are firm practical evidences that Aedes aegypti in Colombia is becoming more and more resistant to insecticides and larvicides (see, e.g.[11,28] and references therein).The latter is also reflected in increasing incidence of registered dengue cases (see Figure 1.1).Therefore, a sole use of chemical substances is not sufficient to control vector population and integrated control strategies can give better results.One of such strategies potentially applicable in Cali, Colombia could be a combination of insecticide control around street catch basins with the use of predacious specious in other mosquito breeding sites, such as ornamental basins filled with stagnant water, located in public areas (shopping and recreation centers, schoolyards, university campuses, parklands, etc.) where the spreading of chemical substances is vetoed by the owners.

Appendix A: Positive Invariance of the Set Ω
In effect, the third equation in (2.1) is logistic; therefore, 0 ≤ x ≤ K x for all initial states with x(0) = x 0 taken from this range.Further, let us consider the projection of Ω onto the plane (m q , m a ).It should be noted that upper bound of m a is obtained according to (3.2).
If we take initial conditions on the left-side vertical edge m q = 0 of twodimensional rectangular region in the space R 2 + = (m q , m a ) illustrated in the Figure A.1, we have that and the associated vectors in the field of directions point rightwards.Alternatively, if the initial conditions are fixed on the right-side vertical edge m q = K q of the rectangular region, it results that meaning that the associated vectors in the field of directions point leftwards.
Further, if we take initial conditions on the lower horizontal edge m a = 0, we have that dm q dt = − (ν + δ q + ϕx) m q < 0, dm a dt = f νm q > 0 and the associated vectors in the field of directions point upwards.Alternatively, if the initial conditions are fixed on upper vertical edge m a = (f νK q ) δ a , it results that dm a dt = f ν (m q − K q ) ≤ 0 while the sign of dm q dt can be either positive or negative.In both cases, the associated vectors in the field of directions point downwards.

Appendix B: Proof of the Theorem 1
Stability analysis is performed using standard technique based on the signs of eigenvalues of the community matrix corresponding to the linearized version of nonlinear system (2.1) 13 which is evaluated at each equilibrium.The community (or Jacobian) matrix of the system (2.1) is given in general terms by .
Further, we proceed to establish the signs of the real parts of all eigenvalues of the Jacobian matrices evaluated at all possible equilibria (3.7), (3.12).

1.
Equilibrium point E 1 defined by (3.7a).At this point, the eigenvalues λ 2 , λ of the Jacobian matrix are the roots of the characteristic polynomial: and one of them is always positive: λ 3 = ρ > 0. The other two eigenvalues, λ 1 and λ (1) 2 , can be either positive or negative.The latter is defined by the sign of free term of the quadratic equation within brackets.If this term is negative (or, equivalently, if Q > 1) then both λ Regarding item (a) it is worthwhile to note that if x(0) = 0 then there exists a two-dimensional stable manifold Ω m = (m q , m a ) ∈ R 2 + \ {(0, 0)} : (B-1) 0 ≤ m q ≤ K q , 0 ≤ m a ≤ f ν δ a K q ⊂ Ω such that for any choice of initial conditions (m q (0), m a (0)) ∈ Ω m the trajectories of the system (2.1) will reach the point E 1 and the total Aedes aegypti population will eventually become extinct 14 .However, it is hardly expected in practice that the mosquito survival threshold without predation, Q, may ever drop below 1.The latter would imply that one adult female mosquito produces, in average, less than one female mosquito and such situation does not comply with the course of nature.
Therefore, item (b) agrees better with common sense and states that E 1 is absolutely unreachable in infinite time (excluding the trivial case when m q (0) = m a (0) = x(0) = 0).

2.
Equilibrium point E 2 defined by (3.7b).The Jacobian matrix evaluated in E 2 can be written using Q (introduced in (3.5)) as . and after some tedious manipulations its characteristic polynomial can be displayed in the following form: Keeping in mind that E 2 exists only if Q > 1, it becomes clear that the above polynomial has one positive root and two roots with negative real part (since the free term of the quadratic equation within brackets is positive).Therefore, E 2 is a saddle point regardless of the value of Q 0 (that is, for either Q 0 < 1 < Q or 1 < Q 0 < Q) and E 2 is conditionally stable.In other words, if x(0) = 0 there exists a two-dimensional stable manifold Ω m (given by (B-1)) such that for any choice of initial conditions (m q (0), m a (0)) ∈ Ω m the trajectories of the system (2.1) will eventually reach the point E 2 ∈ Ω m .The position of E 2 within Ω m is defined by the magnitude of Q > 1.For smaller values of Q > 1 (that is, when mosquitoes barely survive) this point is located closer to the origin and tends to E 1 when Q → 1 + .Alternatively, larger values of Q > 1 will move E 2 towards the right upper vertex of the rectangle Ω m (see Figure A.1).
All roots of this polynomial will have negative real part, if and only if, a free term of the quadratic equation within brackets is positive.The latter implies, in virtue of (3.9), that under the condition Q 0 < 1, the equilibrium point E 3 is locally asymptotically stable.In common terms, it means that if one female mosquito produces, in average, less than one female mosquito due to timely predation of the immature stages, then it is expected that total population of Aedes aegypti will become locally extinct.
Otherwise, if Q 0 > 1, it is concluded that E 3 is a saddle point and there exist a one-dimensional stable manifold such that for (m q (0), m a (0)) = (0, 0) and x(0) ∈ Ω x the trajectories of the system (2.1) will eventually reach the point E 3 .In other words, in absence of immature and adult mosquitoes, the predator will simply reach its carrying capacity using alternative food resources.
We should keep in mind that E 4 exists only if Q 0 > 1.Under this condition, all three eigenvalues given above have negative real parts.Therefore, the equilibrium point E 4 is locally asymptotically stable.

Figure 1 . 1 :
Figure 1.1: Registered cases of dengue reported to Municipal Public Health Secretariat of Cali, Colombia during 2003-2013.
Figure 4.1: Absolute values of sensitivity indexes Υ ν , Υ δq and Υ ϕ = Υ Kx as functions of the product ϕK x for warm-season average temperatures (T =23.9 o C) of Cali, Colombia (see Table 5.1).

Figure 5 . 2 :
Figure 5.2: (a) Level set of Q 0 (ϕ, K x ) for T =23.9 o C; (b) Q 0 (ϕ, K x ) as a real function of the product ϕK x .

Figure 6 . 2 :
Figure 6.2: Population dynamics of immature stages m q (t) during cool, warm, and hot seasons: (a) interaction with an inert predator (Q 0 > 1 for all seasonal temperature variations); (b) interaction with a relatively active predator (Q 0 < 1 for cool and warm seasons, while Q 0 > 1 for hot season).

Figure 6 . 3 :
Figure 6.3: Population dynamics of adult female mosquitoes m a (t) during cool, warm, and hot seasons: (a) interaction with an inert predator (Q 0 > 1 for all seasonal temperature variations); (b) interaction with a relatively active predator (Q 0 < 1 for cool and warm seasons, while Q 0 > 1 for hot season).

Figure A. 1 :
Figure A.1: Vector field directions within projection of Ω onto the plane (m q , m a ) ; the shadowed area also represents the stable manifold Ω m .

2
real parts.If this term is positive (or, equivalently, if Q < 1) then both λ have negative real part.Therefore, (a) if Q < 1 then E 1 is a saddle point; (b) if Q > 1 then E 1 is a nodal repeller.

Table 5 .
1: Parameters of the predator-prey model (2.1) estimated for average temperatures of Cali, Colombia.