A stage structured model of delay differential equations for Aedes mosquito population suppression

Tremendous efforts have been devoted to the development and analysis of mathematical models to assess the efficacy of the endosymbiotic bacterium Wolbachia in the control of infectious diseases such as dengue and Zika, and their transmission vector Aedes mosquitoes. However, the larval stage has not been included in most models, which causes an inconvenience in testing directly the density restriction on population growth. In this work, we introduce a system of delay differential equations, including both the adult and larval stages of wild mosquitoes, interfered by Wolbachia infected males that can cause complete female sterility. We clarify its global dynamics rather completely by using delicate analyses, including a construction of Liapunov-type functions, and determine the threshold level \begin{document}$ R_0 $\end{document} of infected male releasing. The wild population is suppressed completely if the releasing level exceeds \begin{document}$ R_0 $\end{document} uniformly. The dynamical complexity revealed by our analysis, such as bi-stability and semi-stability, is further exhibited through numerical examples. Our model generates a temporal profile that captures several critical features of Aedes albopictus population in Guangzhou from 2011 to 2016. Our estimate for optimal mosquito control suggests that the most cost-efficient releasing should be started no less than 7 weeks before the peak dengue season.


1.
Introduction. Mosquito-borne diseases such as dengue, malaria, and Zika have grown dramatically in recent decades, with over 390 million people infected annually in the tropical and subtropical regions in the world [3,30]. It is largely due to the aggressive spread of the mosquito vectors, combined with increasing urbanization and international travels. As a competent vector for dengue, chikungunya and Zika, Asian tiger mosquito Aedes albopictus has colonized all continents except Antarctica during the last four decades [3,30]. It is the sole vector of dengue in Southern China. In 2014, an unprecedented outbreak of dengue fever occurred in Guangzhou, with 37,354 laboratory confirmed cases [19]. To control dengue, tremendous efforts and funds have been invested in the vector control by larval source reduction and insecticide spraying. However, the ubiquitous emergence of larval sources and insecticide resistance have made these traditional control methods unsustainable [2,38].
The endosymbiotic bacterium Wolbachia has driven the development of a new method for Aedes mosquitoes and disease control since the pioneer study of Xi at al. [31] in 2005. The infection of some Wolbachia strains by Aedes mosquitoes not only restrains the mosquitoes' ability of transmitting disease [28], but also causes the cytoplasmic incompatibility (CI) that wild females crossed with infected males become sterile [8,31]. These critical properties have propelled the development of two parallel disease control strategies: In the first approach called population replacement, Wolbachia infected mosquitoes, male and female mixed, are released in natural areas for a complete invasion of Wolbachia in the wild mosquito population [8,28]. In the second approach called population suppression, only infected males are released in natural areas to suppress the wild population. Since 2015, the largest factory of producing Wolbachia infected male Aedes albopictus in the world has been built up in Guangzhou, where over 97% of Aedes albopictus populations have been suppressed in the treated areas [29].
A deep understanding of mosquito population dynamics can help identify the most influential factors for designing the most effective mosquito control policy. In past few years, we have studied the mosquito population dynamics interfered by Wolbachia, focusing on various aspects such as the population replacement by using delay differential equations [12,13,32,36], the inhomogeneous spatial distribution through reaction-diffusion equations [14,15], the impact of random climate changes through stochastic differential equations [9,10,11], and the maternal transmission leakage of Wolbachia [33,34,37]. However, we included only the terrestrial (adult) stage in most of our models. The three aquatic stages (egg, larva and pupa) were neglected until our very recent studies in an age-stage structure model by discrete equations [35,38]. It has been observed that the density induced intra-specific competition occurs mostly in the larval stage [1,23,26,27]. The competition has found to be a major determinant for the mosquito population growth that elevates mortality rates, delays development time, and influences the female size and fecundity [6,16,18,22,39]. To analyze more thoroughly how the competition affects the efficacy of population suppression, we incorporate the larval stage explicitly in our mathematical model in the current study.
We consider the population dynamics of Aedes albopictus in an isolated area where the mosquitoes are uniformly distributed spatially and evenly distributed in sex. We let A(t) denote the population size of wild adult mosquitoes, and L(t) the population size of larvae. The wild population is interfered by a total number of R(t) male mosquitoes infected by the Wolbachia strains that can induce complete CI in wild females. As usual, we assume first order stage transitions and a first order natural death in both larvae and adults, with the following rate constants: β > 0 -the average number of larvae produced by one female adult from compatible mating with wild males; µ ∈ [0, 1] -the pupation rate of larvae; m > 0 -the natural mortality rate of larvae; α ∈ [0, 1] -the eclosion rate of pupae; δ > 0 -the natural mortality rate of adults. In addition, we let τ 1 > 0 be the average development period from the eclosion of adult to the hatching of first instar larvae of the next generation, and τ 2 > 0 be the average development duration of larvae and pupae.
The number of larvae is increased by the production of females crossed with wild males, and is decreased by death and pupation. For a single female, the probability of a compatible mating with a wild male equals the number of wild males, A/2, divided by the total number of all males, A/2 + R. By combining with the female number A/2 and the rate constant β, we find the production rate of larvae at We follow the idea of the classical logistic model to describe the mortality of larvae [7]: For small L the population decays linearly with the rate m > 0; for large L the larvae compete with each other for the limited resources in the breeding sites and the decay is dominated by a second order term. Let K L > 0 be a carrying capacity type constant characterizing the intensity of the density induced competition. We introduce as the mortality rate of larvae. In our simulation, we assume that K L scales with the size of the mosquito inhabiting area. By combining (1), (2), and the first order pupation, eclosion, and adult death, and integrating the waiting durations τ 1 and τ 2 , we derive the following system of delay differential equations: By using the change of variables supplemented with the conversions of parameters we transform (3) into We will study the dynamics of (6) under the initial conditions for some fixed time t 0 and continuous functions φ(t) and ψ(t).
We give a rather complete description of the global dynamics of (6)- (7) in Section 2 by constructing Liapunov-type functions or using delicate analyses based on the fluctuation lemma. The results are summarized in Theorems 2.1-2.5, and are further interpreted in terms of the original system parameters in Section 3 by reversing (4). Theorem 2.1 shows that the wild mosquito population would not extinct naturally if and only ifᾱβµ > 2δ(m + µ), which defines a threshold relation between the reproduction and mortality of mosquitoes that determines the population's survivorship. Under this condition, Theorem 2.2 shows that if the population is not intervened by Wolbachia. The constants L * and A * can be regarded as the larval and adult carrying capacities for the wild population in the studying area. They both scale with the mosquito inhabiting area, increase in the production rate constants α, β, and µ, and decrease in the mortality rate constants m and δ by nonlinear dependence. Theorems 2.4 and 2.5 give the threshold Wolbachia infected male releasing level A complete population suppression is ascertained when the releasing level R(s) exceeds R 0 uniformly, and is conditional on the initial population size otherwise. When R(s) ≡ R 0 , E 0 (0, 0) is locally stable and there emerges a semi-stable positive equilibrium point E * 2 . When R(s) takes a constant value in (0, R 0 ), E * 2 bifurcates into an unstable E 1 and a stable E 2 for a bistable dynamics of (3). As the mosquito population growth is very sensitive to temperature, we use numerical examples to examine how the temporal profiles of mosquito abundances respond to the variation of the climatic conditions by changing the constant rates of (3) to temperaturedependent rates on a daily basis. Our simulated time-course curves from 2011 to 2016 capture several critical features of natural Aedes albopictus population in Guangzhou, that reproduces quickly from early spring, peaks in later September and early October that coincides with the highest risk season of dengue fever, declines from November and nearly vanishes from December to February. Finally, we use our model to search for an optimal infected male releasing policy to reduce more than 95% of Aedes albopictus population during the high risk season of dengue fever. By analyzing the variation of the least releasing amounts among various control periods, we recommend that the releasing should be started no less than 7 weeks before the target date for the most cost-efficient control.
2. Global stability of the complete suppression state. We study the global dynamics of (6)-(7) in this section. We first claim that the solution (x(t), y(t)) remains positive and bounded in (t 0 , ∞). Indeed, if the positivity fails, then there is t 1 > t 0 or t 2 > t 0 such that x, y > 0 in [t 0 , t 1 ) and x(t 1 ) = 0, or x, y > 0 in [t 0 , t 2 ) and y(t 2 ) = 0. If the first case occurs, then x (t 1 ) ≤ 0 and the first equation of (6) gives and a contradiction. If the second case occurs, then y (t 2 ) ≤ 0 and the second equation of (6) gives y (t 2 ) > 0 and again a contradiction. The boundedness can be proved by a similar method as in the proof of Lemma 2.1 in [13].

2.1.
A condition for natural extinction. Clearly, E 0 (0, 0) is an equilibrium point of (6) and is the only one when R(t) changes in time. One of our major goals is to identify the releasing intensity such that E 0 (0, 0) becomes globally asymptotically stable and the wild mosquito population can be suppressed completely. Interestingly, the following theorem provides a simple theoretical condition under which the wild population would go to extinction naturally. (7) is globally asymptotically stable.
Proof. Let (x(t), y(t)) be a solution of (6), and define the Liapunov-type function From (6) we derive and therefore It follows from b * ≤ 0 that dV /dt < 0, and V (t) decays to zero as t → ∞. Hence (x(t), y(t)) → E 0 as t → ∞, and E 0 is globally asymptotically stable.
Recall from the definition of b and α in (5) we find that The three constants β, µ, and α all measure the reproduction efficiency of female mosquitoes, while the product 2(m + µ)δ in the right side increases in the mortality rates m and δ. Hence (9) gives a threshold relation between the reproduction and mortality of mosquitoes that determines the population's fate between persistence and extinction.

2.2.
The classification of the equilibria. In view of Theorem 2.1, we maintain the following condition in the rest of our discussion without further mention. In the study of the global dynamics of (6), it is useful to classify the equilibria that may emerge when R(t) is identically a constant. It includes, in particular, the important case R(t) ≡ 0 that the wild area is not intervened by Wolbachia infected males. Let R(t) ≡ R ≥ 0 be a constant. If E(x, y) is an equilibrium point of (6), then y = αx and x(1 + x) = by 2 /(y + R). It follows that

MUGEN HUANG, MOXUN TANG, JIANSHE YU AND BO ZHENG
The discriminant of g(x, R) is If we define then We use these calculations to determine the existence of equilibria with positive coordinates (or positive equilibria in short) of (6) besides the universal equilibrium point E 0 (0, 0).
Both x 1 and x 2 , whose expressions follow from (11) directly, are well-defined since ∆ R > 0 in this case. Also, since we have αb * − R > 0 and so x 2 > 0. As g(0, R) = R > 0, it is also seen that x 1 > 0.

2.3.
The global dynamics of the wild mosquito population. When R(t) ≡ 0, (6) describes the dynamics of the wild mosquito population without the intervention of Wolbachia. The discussion in Section 2.2 indicates that (6) has two nonnegative equilibria: E 0 and E * 1 (b * , αb * ). We show that E * 1 is globally asymptotically stable. Thus b * and αb * define the respective carrying capacities of larvae and adults of the wild population in the studying area.
Proof. We only need to prove the global asymptotical stability of E * 1 . Let (x(t), y(t)) be a solution of (6)- (7). We first prove In fact, if (17) is not true, then either x = 0 or y = 0 by the positivity of x and y. We claim that x = 0 and y = 0 could only occur simultaneously. Suppose y = 0. Then there is an infinite sequence {t n } such that y(t n ) → 0 and y (t n ) ≤ 0. Substituting t = t n in the second equation of (6) yields αx(t n − τ 2 ) ≤ y(t n ), and so x = 0. Similarly, when x = 0, there is an infinite sequence {s n } such that x(s n ) → 0 and x (s n ) ≤ 0. Taking limit in the first equation of (6) along {s n } gives and so y = 0. This verifies the claim.
To complete the proof of (17), we proceed with the assumption that x = y = 0. We construct a sequence, denoted again by {s n }, in the following way: Let s 1 > t 0 be the least time t at which x(t) equals a half of the minimum value of φ(t) over In general, we take s n as the least time in (s n−1 , ∞), n ≥ 2, such that By the construction, it is seen that s n > s n−1 , Inserting the first part of (18) into the first equation of (6) gives and so y(s n − τ 1 ) ≤ x(s n )(1 + x(s n ))/b. At t = s n − τ 1 , the second equation of (6) gives for sufficiently large n. This shows that y(t) increases strictly at each time t = s n − τ 1 . Since y = 0, there must be a sequence of positive integers n 1 < n 2 < · · · such that y(t) decreases in a subinterval within ( Evaluating y (t i ) by the second equation of (6) gives As s ni > s ni − τ 1 > t i > t i − τ 2 , the second part in (18) leads to

MUGEN HUANG, MOXUN TANG, JIANSHE YU AND BO ZHENG
Now, evaluating x (s ni ) by the first equation of (6) gives Provided that i is sufficiently large, this would imply x (s ni ) > 0, which contradicts (18) and completes the proof of (17).
Having shown (17), the rest of the proof follows a standard approach. In our discussion below, we will use increasing sequences {s n } and {t n } frequently that may differ in different contents with {s n } → ∞ and {t n } → ∞ as n → ∞. By the fluctuation lemma, see Lemma A.1 in [24], we may find {s n } and {t n } such that x(s n ) → x, x (s n ) → 0, and y(t n ) → y, y (t n ) → 0, as n → ∞. Taking the limit in the second equation of (6) along {t n } yields Similarly, taking the limit of the first equation of (6) along s n gives Then x ≥ b * , and y ≥ αb * . On the other hand, let {s n } and {t n } be such that x (s n ) → 0, and y (t n ) → 0 as n → ∞. Then the same limiting processes above yieldȳ ≤ αx andx(1 +x) ≤ bȳ ≤ αbx, leading tox ≤ b * andȳ ≤ αb * . Hence x =x = b * , and y =ȳ = αb * , that confirm the global asymptotical stability of E * 1 (b * , αb * ).
2.4. The bistable dynamics under insufficient control. When the Wolbachia infected males are released at a low level that R(t) ∈ (0, R * ) (see (12) for the definition of R * ), we show that E 0 is locally, but not globally, stable. It indicates that a complete population suppression is conditional on the initial population size.
For the case R(t) ≡ R ∈ (0, R * ), the discussion in Section 2.2 shows that (6) has three nonnegative equilibria E 0 , E 1 (x 1 , y 1 ), and E 2 (x 2 , y 2 ) with x 1 and x 2 being defined in (14). In this case, our next theorem reveals that system (6) displays a bistable dynamics with two stable equilibria E 0 and E 2 .
2.5. The threshold releasing level. We finally show that R * defines a threshold releasing level of Wolbachia infected male mosquitoes for a complete suppression of wild populations: When inf [t0−τ,∞) R(t) > R * , both larval and adult populations go to extinction; when R(t) ≡ R * , the positive equilibrium point E * 2 (x * 2 , y * 2 ), where x * 2 = √ αb − 1 and y * 2 = αx * 2 , becomes semi-stable. Theorem 2.5. Let (x(t), y(t)) be the solution of (6)-(7). Then we have the following statements. ( Proof. (1) We omit the proof since it is virtually identical to that of Theorem 2.4.
3. Application to population suppression.

3.1.
Interpreting the theorems in original system parameters. By reversing the change of variables defined in (4), we have It can convert (6)  We interpret L * and A * as the larval and adult carrying capacities for the wild population in the studying area. They both scale with the carrying capacity constant K L , increase in the production rate constants α, β, and µ, and decrease in the mortality rate constants m and δ by nonlinear dependence. We note that a higher pupation rate µ of larvae induces a larger reduction of larval numbers temporally, but the larger reduction is compensated by the production of more larvae through the increase of adults at steady-state. The threshold releasing level R * defined in (12) is converted to Theorem 2.5 shows that if R(s) exceeds R 0 uniformly, then the wild population will be wiped out successfully. Otherwise, the detailed analyses summarized in Theorem 2.4 indicate that a complete population suppression is conditional on the initial population size. As K L scales with the size of the mosquito inhabiting area, (24) reveals that the threshold releasing level R 0 also scales with the inhibiting area as we conceive intuitively. Like L * and A * , R 0 also increases in α, β, and µ, and decreases in m and δ. We denote again by E with the same subscript or superscript the equilibria of (3) as those for (6). Clearly, E 0 (0, 0) is also a universal equilibrium point of (3) and is the only one when R(s) is not a constant. When R(s) is identically a constant R, the existence of additional equilibrium points can be classified as in Section 2.2 by converting the coordinates (x, y) into (L, A) according to (20). When R ∈ (0, R 0 ), (3) has two positive equilibria E 1 (L 1 , A 1 ) and E 2 (L 2 , A 2 ) corresponding to E 1 (x 1 , y 1 ) and E 2 (x 2 , y 2 ) defined in (14). Theorem 2.4 shows that E 0 and E 2 are locally stable while E 1 is unstable, and (3) displays a bistable dynamics. When R = R 0 , (3) has a unique positive equilibrium point E * 2 (L * 2 , A * 2 ) corresponding to E * 2 (x * 2 , y * 2 ) defined by (16). Theorem 2.5 shows that E * 2 is stable from right side and unstable from left side, and so it is also called semi-stable.
To demonstrate more specifically how our mathematical results can be applied in the control of Aedes albopictus in Guangzhou, we fix the following parameter values β = 2, m = 0.07, µ = 0.1,ᾱ = 0.95, δ = 0.1, τ E = 4, τ L = 6, τ P = 2, τ A = 16, (25) in accordance with the life table of Aedes albopictus in Table 1. From (25), we have the estimatesτ 1 ≈ τ E + τ A /2 = 12, andτ 2 = τ L + τ P = 8. In the simulations, we take K L = 5 × 10 5 as an example, which can be replaced by any larger integer without changing the system dynamics. It gives the threshold releasing level R 0 ≈ 1.073 × 10 6 . Take for example R = 10 5 ∈ (0, R 0 ). Then (3)  By Theorem 2.4, E 1 is unstable, while E 0 and E 2 can attract solutions with the initial data (φ(s), ψ(s)) lying within the appropriate intervals. In Figure 1A, a solution curve (L, A) is sketched for (φ, ψ) = (1.1×10 5 , 7.5×10 3 ), a case not covered by Theorem 2.4 since φ > L 1 and ψ < A 1 . The solution displays an interesting and complex dynamical behavior that rotates around the unstable equilibrium point E 1 in several circles before moving towards E 0 . This example also shows that E 0 may attract solutions with the initial data not covered by Theorem 2.4, which is further evidenced by the solution with (φ, ψ) = (10 4 , 6 × 10 4 ) in Figure 1B as φ < L 1 and ψ > A 1 . Similarly, as shown in Figure 1C, E 2 can also attract solutions with φ > L 1 and ψ < A 1 such as (φ, ψ) = (2 × 10 5 , 10 4 ), or solutions with φ < L 1 and ψ > A 1 such as (φ, ψ) = (10 4 , 2 × 10 5 ). When R = R 0 = 1.073 × 10 6 , E 1 and E 2 coincide to E * 2 (L * 2 , A * 2 ) with L * 2 = 1.6562 × 10 6 , A * 2 = 1.5734 × 10 6 . As shown in Figure 1D, similar observations can be made about the attractivity of E 0 and E * 2 : The solutions with φ > L * 2 and ψ < A * 2 such as (φ, ψ) = (3×10 6 , 5×10 5 ) or (10 7 , 10 6 ), and the solutions with φ < L * 2 and ψ > A * 2 such as (φ, ψ) = (5 × 10 5 , 1.6 × 10 6 ) or (10 6 , 3 × 10 6 ), can approach either E 0 or E * 2 ultimately. It remains a rather subtle question to determine the ultimate fate of (L, A) when φ and ψ take constant values that are not covered by our Theorems 2.4 and 2.5. It is more challenging to discuss the general case when φ and ψ take non-constant values that cross over (L 1 , A 1 ) or (L * 2 , A * 2 ). 3.2. The impact of temperature on population suppression. The mosquito population growth is very sensitive to climatic conditions such as temperature and rainfalls. Guangzhou has a typical subtropical climate with a hot and rainy summer and a warm and dry winter. July and August are the hottest months with the temperature higher than 30 o C in most days, and January is normally the coldest month but the temperature is rarely less than 10 o C (http://www.cma.gov.cn/). The mosquito density in Guangzhou peaks in late September and early October, which overlaps with the high-risk period of dengue fever [18,21,22]. Either E0 or E2 may attract solutions with φ > L1 and ψ < A1, or φ < L1 and ψ > A1; see (26) for the values of L1 and A1. D. Either E0 or E * 2 may attract solutions with φ > L * 2 and ψ < A * 2 , or solutions with φ < L * 2 and ψ > A * 2 .

Para. Value
Reference µ E ρ E = 0.24, a E = 10798, b E = 100000, c E = 14184 [6,17,25] µ = µ L ρ L = 0.2088, a L = 26018, b L = 55990, c L = 304.6 [6,17,25] α = µ P ρ P = 0.384, a P = 14931, b P = −472379, c P = 148 [6,17,25]   Understanding the vector-environment relationship is essential for the control of mosquito populations and the prevention of mosquito-borne diseases. To integrate temperature into our model, we change the constant mortality rates m and δ in (3) to the rates depending on the temperature T as specified in Table 2. In addition, we modify the production rate constants in (3) by introducing the temperaturedependent developmental rates of egg (µ E ), larva (µ L ), and pupa (µ P ) according to the following formula given in [6]: , for S = E, L, P, In (27), K is the mean temperature in Kelvin, and γ = 1.987 cal K −1 mol −1 is the universal gas constant. ρ S is the development rate at 298 Kelvin assuming no enzyme inactivation, and c S is the temperature when half of the enzyme is deactivated. Their values, along with a S and b S , are listed in Table 2.  (3) with R ≡ 0, supplemented by the temperature-dependent rates estimated from (27) and Table 2.
By using the temperature data from 2011 to 2016 in Guangzhou, we can use (27) and Table 2 to estimate the temperature-dependent mosquito production and mortality rates on a daily basis, and then use (3) with R ≡ 0 to simulate the daily Aedes albopictus abundances in these years. We set January 1, 2011 as the initial time t 0 = 0, and fix the initial number of larvae at 10 6 and the number of adults at 10 4 over the appropriate initial time period, corresponding to the constant K L = 5×10 5 as in Figure 1. As shown in Figure 2, the simulated time-course curves exhibit a strong periodicity as the temperature in Guangzhou, and capture several critical features of natural Aedes albopictus population in Guangzhou [18,21,22]. The wild mosquito population reproduces quickly from early spring, peaks in late September and early October, declines from November and nearly vanishes from December to February in next year. As the system parameters change in time, the threshold releasing level R 0 defined in (24) may vary accordingly from day to day. As shown in Figure 3, R 0 exhibits a quasi-periodicity as temperature annually, and takes large values from July to October, that coincides with the peak season of mosquito abundances and disease transmissions. Hence the period from July to October is the critical control season of Aedes albopictus population and dengue fever.  (24) with the daily rates estimated from (27) and Table 2. It exhibits a quasi-periodicity as temperature annually, and peaks from July to October during the high-risk season of dengue fever.  Table 3. The least releasing levels r1 in phase I, r2 in phase II, and the total numbers in the two phase to reduce > 95% of wild Aedes albopictus on September 12, 2012, and keep the same low level in phase II of the five weeks behind. The same parameter values and the initial data in Figure 2, except R(t) = r1 in phase I, and R(t) = r2 in phase II, were inserted into (3) for simulations.

3.3.
Mosquito suppression during the peak season. As shown in Figure 2, the highest Aedes albopictus density in 2011 and 2012 in Guangzhou was observed around September 12, which was also at the beginning of the highest risk season of dengue fever. An efficient strategy to control dengue fever is to reduce the wild mosquito population to a substantially lower level within the dengue high-risk season. We use our model to help design a Wolbachia infected male releasing policy, targeting at a more than 95% reduction on September 12, and keeping the same adult mosquito density from September 12 to October 18 for 5 weeks. We divide the releasing into two phases, a strong suppression phase with a large releasing amount (phase I) before September 12, followed by the maintenance phase in the five weeks after September 12 with a small releasing amount (phase II). We assume that the infected males are released once in three days at the same amount r 1 in phase I and 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 Figure 4. Suppression of the wild Aedes albopictus population during the high-incidence season of dengue fever in 2012. The curves were generated by substituting the same parameter values and the initial data described in Table 2 into (3). The mosquito population is reduced more than > 95% on September 12, 2012, and the same low level is kept in phase II of the five weeks behind. r 2 in phase 2. Although the total number R(t) may change from day to day, the number of infected males actively competing with wild males for mating may not change considerably. Thus we assume R(t) ≡ r 1 (or r 2 ) within the two phases in our simulation. We test the releasing policy by setting the starting dates in 5 to 10 weeks before September 12 and search for the most cost-efficient starting date. We set the target date on September 12, 2012 to demonstrate the test. We use the same set of parameter values and the initial data as in the simulation of Figure  2, with the only exception of changing R from zero to r 1 in phase I from the starting date to the target date, and to r 2 in the second phase. By numerical simulations, we estimate the least values of r 1 and r 2 such that the mosquito abundance is reduced to the target level mentioned above. The values corresponding to the starting dates in 5 to 10 weeks before the target date are listed in Table 3. The temporal profiles of the adult daily abundance in 2011 and 2012 are shown in Figure 4. It is seen from Table 3 that r 1 decreases sharply from 4.4 × 10 9 for the five week releasing to 5.11 × 10 8 (less than 1/8 of 4.4 × 10 9 ) for the six week releasing in phase I. It is further reduced from 5.11 × 10 8 to 1.48 × 10 8 for the seven week phase I. The reduction of r 1 is less significant when the length of phase I is shifted to 8, 9, and 10 weeks. The change of r 2 in different phase I periods is insignificant as we expect. However, the total releasing numbers show about the same reduction pattern as r 1 . By analyzing the variation of the least releasing amount r 1 and the total numbers among different phase I durations, we recommend the mosquito control to be started no less than 7 weeks before the target date.