Mathematical Biosciences 258 (2014) 77–84

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Review

Multiobjective Genetic Algorithm applied to dengue control Helenice O. Florentino∗, Daniela R. Cantane, Fernando L.P. Santos, Bettina F. Bannwart UNESP – Univ. Estadual Paulista, IBB – Instituto de Biociências, Departamento de Bioestatística, Distrito de Rubião Júnior, s/n, Botucatu, SP, Brazil

a r t i c l e

i n f o

Article history: Received 29 November 2013 Revised 14 August 2014 Accepted 22 August 2014 Available online xxx Keywords: Dengue Multiobjective optimization Genetic algorithm

a b s t r a c t Dengue fever is an infectious disease caused by a virus of the Flaviridae family and transmitted to the person by a mosquito of the genus Aedes aegypti. This disease has been a global public health problem because a single mosquito can infect up to 300 people and between 50 and 100 million people are infected annually on all continents. Thus, dengue fever is currently a subject of research, whether in the search for vaccines and treatments for the disease or efficient and economical forms of mosquito control. The current study aims to study techniques of multiobjective optimization to assist in solving problems involving the control of the mosquito that transmits dengue fever. The population dynamics of the mosquito is studied in order to understand the epidemic phenomenon and suggest strategies of multiobjective programming for mosquito control. A Multiobjective Genetic Algorithm (MGA_DENGUE) is proposed to solve the optimization model treated here and we discuss the computational results obtained from the application of this technique. © 2014 Elsevier Inc. All rights reserved.

1. Introduction The dengue virus is most commonly transmitted by the mosquito vector Aedes aegypti, but can be transmitted by other members of the genus Aedes including Aedes albopictus. The tropical regions are most affected by dengue fever due to environmental conditions that favor the development and proliferation of A. aegypti. In hot and rainy periods there occurs increased density of mosquitoes, causing a greater incidence of dengue fever. This disease is caused by four distinct serotypes known as DEN1, DEN2, DEN3 and DEN4. If a person is infected by one of the four serotypes, he will never again be infected by the same serotype (homologous immunity), but loses immunity to the three other serotypes (heterologous immunity) in about 12 weeks and then becomes more susceptible to developing dengue hemorrhagic fever [4,6,16]. A female mosquito becomes infected with dengue by feeding on the blood of an infected person and then transmits the virus to other people through subsequent blood feedings. A mosquito remains infected for the duration of its life and females can also transmit the virus to their offspring. According to Health Surveillance Secretariat of Brazil, in 45 days of life, a single mosquito can infect up to 300 people. The incubation of the virus in humans ranges from 3 to 15 days after the bite of the vector, with an average of 6 days [4,16,21].

∗

Corresponding author. Tel.: +55 14 3880 0070. E-mail address: [email protected] (H.O. Florentino)

http://dx.doi.org/10.1016/j.mbs.2014.08.013 0025-5564/© 2014 Elsevier Inc. All rights reserved.

There is an estimate by the World Health Organization (WHO) that annually from 50 to 100 million people are infected in more than 100 countries on all continents. With nearly half of the world’s population at risk of dengue infection, the lack of specific treatment or vaccine increases the severity and risk. It is estimated that 500,000 people with dengue hemorrhage fever (DHF) require hospitalization each year – a large proportion of these are children. It is fatal in approximately 2.5% of those affected [21]. The prevention and control of dengue fever depend on measures that effectively control the mosquito A. aegypti. In Brazil, the main controls currently used by the Superintendence of Endemic Disease Control (SUCEN, Superintendência de Controle de Endemias) are: physical control, done by public health officials and residents through eliminating breeding; and chemical control, done by applying insecticides. Insecticides play an important part in controlling dengue vectors. However, over the time the vectors acquire resistance, requiring an increase of chemicals, increasing the cost of control and affecting public health [20]. It is very difficult to control or eliminate A. aegypti mosquitoes with these control measures. This is due to the fact that they have adaptations to the environment that make them highly resilient. The mosquito has the ability to bounce back rapidly to initial numbers after disturbances resulting from natural phenomena (e.g., droughts) or human interventions. Thus, the control of the mosquito that transmits dengue is the subject of much research. There are several works with the proposal of assisting in the search for optimal control of the mosquito that transmits dengue. Caetano and Yoneyama [5] present an application of the optimal control theory in dengue epidemics. In this work, the dynamics of this

78

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

insect-borne disease are modeled as a set of non-linear ordinary differential equations. These include the effect of educational campaigns organized to motivate the population to break the reproductive cycle of the mosquitoes by avoiding the accumulation of standing water in open-air receptacles. The functional cost is such that it reflects a compromise between the health of the population and the actual spending on insecticides and educational campaigns. Rodrigues et al. [18] present a dynamic model which describes the dynamics of the dengue mosquito, the number of infected individuals and people’s motivation to combat the mosquito. The functional cost depends on the costs of medical treatment of the infected people as well as on the costs related to education and sanitation campaigns. They present two approaches for solving the problem: one using the optimal control theory and the other carried out by first defining the problem and then solving it with nonlinear programming. Other authors have proposed approaches of optimal control which are applied to dengue differently than the approach that will be presented in the present work ([1,2,20] and others). Brazil has started biological control with genetically modified mosquitoes to fight dengue. In this control, male mosquitoes that have been rendered sterile by irradiation in laboratory are used. This technique is called the Sterile Insect Technique (SIT) and has been shown to be effective in controlling agricultural pests, such as the blowfly eradication program. The sterile males compete with the wild males for female insects. Adult females mate only once in their lifetime and, after mating, reject other males. Additionally, the female lays about 100 eggs. The release is performed gradually and in different places to increase the probability of the birth of mosquitoes. However, if a female mates with a sterile male it will have no offspring, thus reducing the next generation’s population [3,8,14,17]. The main aim of the present work is to study techniques of multiobjective optimization that assist in solving problems involving the control of the mosquito that transmits dengue fever. We do not include all biological and mathematical details in order to give greater importance to the multiobjective technique proposed. The population dynamics of the mosquito that transmits dengue fever is studied in order to understand the epidemic phenomenon and suggest new multiobjective programming strategies for the control of the mosquito. A Multiobjective Genetic Algorithm (MGA_DENGUE) is proposed to solve the optimization model treated here and we discuss the computational results obtained from the application of this technique. 2. Mathematical model Esteva and Yang [9] proposed a mathematical model for analyzing the effect of the introduction of sterile male insects in the environment for control of the mosquito that transmits dengue. In this proposal, the mosquito’s life cycle was considered as being divided in two phases: the aquatic phase, when mosquitoes are in the egg, larva and pupa stages; and the terrestrial phase, when mosquitoes are at the adult stage. The size of the mosquito population in the aquatic phase at time t is denoted by A(t). Mosquitoes in the aquatic phase enter the terrestrial phase at a γ per capita rate, with r proportion of the females and (1 − r) of the males. In the terrestrial phase, the wild adult mosquitoes were divided into three compartments: immature females (before mating), I; natural (or wild) male insects, M; and fertilized females, F. By inserting sterile mosquitoes into the system, it creates two more compartments: sterile (irradiated or transgenic) insects, S; and unfertilized females (removed female), U, which are the females that mate with sterile males. The flows from immature females I to fertilized females F and unfertilized females U compartments depend mainly on the number of encounters with natural and sterile males. Then the probability of an encounter of a female I with natural male insects is given by

(M/(M + S)) and the per capita rate at which female insects are fertilized is (β M/(M + S)), where β is the mating rate of the insects. The probability of an encounter of a sterile male S with a female I is given by (pS/(M + S)), where p, 0 ࣘ p ࣘ 1, is the proportion of sterile insects that are released in the appropriate places. The effective fertilization during mating could be diminished due to sterilization, which leads one to assume that the effective mating rate of sterile insects is given by qβ , 0 ࣘ q ࣘ 1. Putting together the assumptions above, we get that (β S S/(M + S)), were β S = pqβ . Finally, α is the rate at which the population of the sterile mosquitoes S is sprayed in the environment. The net oviposition rate per female insect is proportional to their density, but it is also regulated by a carrying capacity effect depending on the occupation of the available breeder sites. It is assumed that the per capita oviposition rate is given by φ (1 − A/C) where φ is the intrinsic oviposition rate and C is the carrying capacity related to the amount of available nutrients and space. The per capita mortality rates of the aquatic phase, immature females, fertilized females, unfertilized females and sterile male insects are denoted by μA , μI , μF , μU , μM and μS respectively. The dynamics of this system are illustrated in Fig. 1. The parameters A(t), S(t), I(t), M(t), U(t) and F(t) represent the size of each segment of the mosquito population at time t: In addition to the biological control with sterile mosquitoes, we also consider chemical control with insecticide, which has played an important role in controlling A. aegypti. The insecticide acts only in the adult stage of the mosquito and not in the aquatic phase. In chemical control, the amount and cost of insecticide used is generally proportional to the area treated and not to the size of the population. However, the damage done by mosquitoes is proportional to the number of individual per unit area. Thus, the cost/benefit ratio of insecticidal control increases as the size of the population increases. In sterile male release programs, the number of sterile male released and their cost is proportional to the size of the population. The cost/benefit ratio decreases with a decrease in size of the population. Then, it can be wise to use integrated insecticidal methods with sterile male technique. Considering the control variables u1 = u1 (t) and u2 = u2 (t), associated respectively with chemical control and biological control, we present the mathematical model (1) for the population dynamics using the insecticide and sterile mosquitoes controls as proposed by Thome et al. [20].

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

  dA A =φ 1− F − (γ + μA )A dt C  βM βS S dI = rγ A − + + μI + u 1 I dt M+S M+S

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

dM = (1 − r)γ A − (μM + u1 )M dt

β MI dF = − (μF + u1 )F dt M+S

(1)

dS = u2 − (μS + u1 )S. dt

The state variable U, related to removed female mosquitoes, is decoupled from the dynamic system described by

dU βS SI = − μU U − u1 U. dt M+S To study the system state (1), Ref. [20] used the initial conditions presented in (2) because these are the equilibrium conditions

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

Aquatic phase

79

Terrestrial phase Fig. 1. Mosquito population dynamics with biological control with sterile males.

of the system and the worst possible situation from the standpoint of controlling dengue. More details can be found in Refs. [10,20].

⎧ C (R − 1) ⎪ ⎪ A(0) = A0 = ⎪ ⎪ R ⎪ ⎪ ⎪ ⎪ ⎪ r γ A ⎪ 0 ⎪ ⎪ I (0) = I0 = ⎪ ⎪ μI + β ⎪ ⎨ (γ + μA )CA0 F (0) = F0 = ⎪ ⎪ ⎪ φ(C − A0 ) ⎪ ⎪ ⎪ ⎪ (1 − r)γ A0 ⎪ ⎪ ⎪ M(0) = M0 = ⎪ ⎪ μM ⎪ ⎪ ⎪ ⎩ S(0) = S0 = 0,

0 ࣘ t ࣘ T. Therefore, the model is:

Minimize J(u1 , u2 ) = [J1 , J2 , . . . , Jk ] (u1 (t),u2 (t))

Subject to System (1). Initial conditions (2). (2)

where R = φ rγ β /(γ + μA )(β + μI )μF define the average number of secondary female insects produced by a single female insect. Esteva and Yang [9,10] propose a full discussion about this ratio, showing that R >1 is a stable endemic equilibrium. 3. Multiobjective optimization problem Thome et al. [20] discussed the high cost of chemical and biological controls, as well as the cost–benefit of each of these. The authors proposed an optimal mono-objective control model to determine u1 and u2 that optimizes the cost involved in this process – see Ref. [20] for more information. The fact is that the optimization problem addressed here is a multiobjective problem that presents many conflicting objectives: to minimize costs of sterile mosquitoes and insecticide; to maintain the maximum number possible of the sterile mosquitoes placed on site; to minimize the number of insects in the aquatic and adult phases (fertilized females and natural males). Thus, multiobjective optimization techniques are appropriate for this problem. We propose the multiobjective model (3) to determine the optimal control u1 and u2 , where the decision variables u1 = u1 (t) and u2 = u2 (t) are related to chemical and biological controls, respectively. Our aim is to minimize the costs involved in mosquito control and decreasing the populations of natural insects, be they economic, social or environmental costs. What we want is to determine u1 and u2 that optimize some costs involved in chemical and biological control processes in period

(3)

Final conditions over some variable system state (1). Non-negativity conditions : u1 (t) ≥ 0 and u2 (t) ≥ 0, all 0 ≤ t ≤ T. The multiobjective optimization model can have from 2 to k objectives, according to the preference of the user and the characteristics of the problem, where k is a real limited integer. Some objectives that can be used in this multiobjective optimization model are: 1. To minimize the investment with insecticide;

Minimize J1 =

T 0

(w1 u1 )dt

where 0 ࣘ u1 (t) ࣘ 1 for all t in the considered time period 0 ࣘ t ࣘ T. w1 is a weight or to the monetary cost associated with insecticide. T We can also define J1 as J1 = 0 (w1 u21 )dt. The square in the variable amplifies the effects of large variations in the minimization process. 2. To minimize the investment in the production of sterile mosquitoes;

Minimize J2 =

T 0

(w2 u2 )dt.

where u2 (t) ࣙ 0 is the control variable related to the insertion of sterile mosquitoes in the environment. w2 is a weight or can be the cost of introducing sterile male mosquitoes in the environment. To amplify the effects of large variations, we can also define J2 = T 2 0 (w2 u2 )dt. 3. To minimize the number of the fertilized females;

Minimize J3 =

T 0

(w3 F )dt.

80

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

This is a social cost because the fertilized females are responsible for carrying and transmitting the dengue virus and for producing new individuals in the population. w3 is a weight or to the monetary cost associated with the treatment of dengue. We can square F to amplify the effects of large variations. 4. To maintain the number of the sterile mosquitoes that remain in the environment.

Minimize J4 = (−1)



T 0

(w4 S)dt.

where the functional J4 is related to the loss of the sterilized mosquitoes due to insecticide. The insertion of sterile mosquitoes is an expensive process, so we hope that the insecticide kills the minimum number possible of the insects inserted in the environment, thus this objective can be rewritten by removing the (−1) from the function and changing “Minimize” to “Maximize”. w4 is a weight or can be the importance given to biological control or the cost of introducing sterile male mosquitoes in the environment. 5. We can also group one or more objectives in one function; for example, to minimize the insecticide cost and to minimize the cost of sterile mosquitoes, or to minimize the amount of fertilized females in the system (social cost) and to preserve the maximum possible number of introduced sterile mosquitoes, etc. We propose a genetic algorithm (MGA_DENGUE) to solve the optimization problem (3). The main advantages of this method is simplicity of computational programming; the fact of not requiring changes in system (1); the ease of adding constraints and solving the problem only once, since the MGA_DENGUE offers several solutions and allows for the evaluation of the entire search space. 4. Multiobjective Genetic Algorithm The concept of genetic algorithms (GA) was developed by Holland and coworkers in the 1960s and 1970s [13]. GA were inspired by the evolutionist theory for explaining the origin of species. The strong individual (an individual that has good fitness) in a population has a greater chance of passing its genes to future generations via reproduction (recombination or crossover). Species carrying the correct combinations in their genes become dominant in their populations. Sometimes, random changes may occur in genes (mutations) and new species evolve from the old ones (evolution). Unsuccessful changes are eliminated by natural selection (selection). In GA terminology, a solution is called individual and made of discrete units. A collection of individuals is called a population. The population is normally randomly initialized. GA use the operators to generate new solutions from existing ones: selection, crossover and mutation. In crossover, two solutions are generally combined to form a new individual. The solutions are selected among existing solutions in the population with a preference for fitness, so that new individual is expected to inherit good characteristics. By iteratively applying the crossover operator, characteristics from good individuals are expected to appear more frequently in the population, eventually leading to convergence to an overall good solution. The mutation operator introduces random changes in the characteristics of individuals, generally applied at the discrete unit of solution level if the probability of changing the properties of a unit is very small – typically less than 1%. GA have been used to solve mono-objective and multiobjective optimization problems. Multiobjective optimization is very different from single-objective optimization because in mono-objective problems one attempts to obtain the best design or decision, which is usually the global optimum. In case multiobjective there does not exist a single solution which is the best with respect to all objectives,

Fig. 2. Illustration of dominance and Pareto front regarding two objectives in an optimization model.

there exists a set of solutions which are superior to the rest of the solutions in the search space when all objectives are considered. These solutions are called Pareto-optimal solutions or efficient solutions. The rest of the solutions are known as dominated solutions. For illustrations, let us consider two feasible solutions: u∗ = ∗∗ (u∗1 , u∗2 ) and u∗∗ = (u∗∗ 1 , u2 ) for the multiobjective problem (3) with two objectives (minimize (J1 ) and maximize (J2 )) and suppose that J1 (u∗ ) < J1 (u∗∗ ) and J2 (u∗ ) > J2 (u∗∗ ) i.e. u∗ is better solution than u∗∗ because u∗ improve all objectives of the problem as compared with u∗∗ . Then, u∗∗ is strictly dominated by u∗ (see Fig. 2). If u∗ is not dominated by any other feasible solution, then u∗ is an efficient solution or a Pareto-optimal solution. The set of all solutions that are not dominated by any other solution is called the Pareto-optimal solutions set or the efficient solutions set. The values associated with all objectives corresponding to the Pareto-optimal solutions set can be represented in objective space and this representation defines the Pareto front. In Fig. 2, the point ( J1 (u∗ ), J2 (u∗ )) belongs to the Pareto front. A solution for the problem can be chosen over the others by a designer (called the decision maker) according to his preference or necessity. The multiobjective GA generates sets of solutions that iteratively allow for computation of an approximation of the entire Pareto front. The first multiobjective GA, called Vector Evaluated Genetic Algorithms (or VEGA), was proposed by Schaffer [19]. Afterward, another important multiobjective evolutionary algorithm, known as the MOGA (Multiobjective Genetic Algorithm), was developed [11]. Currently, there are many variations of multiobjective GA in the literature. Several are well-known and credible algorithms and have been used in many applications, and their performances were tested in several comparative studies [7]. We propose Multiobjective Genetic Algorithm (MGA_DENGUE) to solve the optimization problem (4). In this algorithm, a population of solutions is generated (individuals) (u1 (t), u2 (t)) for the problem (3) and this population evolves according to genetic operators (selection, crossover and mutation) in order to promote a tendency for solutions to represent efficient or non-dominated solutions which are increasingly effective and more varied as the evolutionary process continues. The solutions of the population are defined by matrices containing 2 lines and f + 1 columns, or T + 1 columns when t = 0, 1, 2, . . . , T, where T is the total period of the application of the controls. The first line refers to the values of the decision variable u1 and the second is related to u2 . Each column represents the discrete period T from initial time t0 to the final time tf . Thus, each element (i, j) in this matrix (individual), i = 1, 2 and j = 0, 1, . . . , f, represents a value for the decision variable in time tj , as illustrated in Fig. 3. The initial population is generated randomly with n individuals, in the form of Fig. 3, and in MGA-Dengue we use 0 ࣘ u1 (t) ࣘ 1 and 0 ࣘ u2 (t) ࣘ 10M(t) for all t in the considered time period, or 0 ࣘ t ࣘ T. We use this range of values for u2 because it is advisable in SIT

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

t0=0

t1=1

t2=2

...

tf=T

u1(t)

u1(0)

u1(1)

u1(2)

...

u1(T)

u2(t)

u2(0)

u2(1)

u2(2)

...

u2(T)

penalized by multiplying the level by a positive integer h = 1000. We added a constraint on the state variable F:

F (t) ≤ Ffix

Fig. 3. Solution structure or individual.

Fig. 4. Level of dominance and fitness for feasible solutions in the objective space considering two objectives in the optimization model.

technique to insert an amount of sterile males 5–10 times the number of natural males to create greater competition. In each iteration, all individuals of the population are evaluated and some individuals are randomly selected to procreate, giving preference to the best one. The worst ones are replaced with new ones generated based on the genetic operators, as described next. The evaluation of each individual of the population is made based on its fitness. The fitness is calculated by creating a ranking originating from the level of individual dominance. The solution belonging to level 1 are determined by finding all nondominated solution for the population (those belonging to the Pareto front). The solutions for level 2 are the best nondominated solutions found when the solutions for level 1 are removed from the population, and so on for all levels (see Fig. 4). The fitness of each individual of the population is calculated by

fi =

1 , h Li

81

(4)

where fi is the fitness of the individual i, Li is the level of the individual i and h is a real number such that h = 1 if the solution i is feasible (if it satisfies all the constraints of problem (3)) and h > 1 in a contrary case. h can be seen as a penalty on unfeasible solutions. We used a large value for h in the case of an unfeasible solution (h = 1000). In this paper, the lower the level, the better the fitness, because the solutions with level equal to 1 are the best solutions. This solution group is an approximation of the Pareto front. Infeasible solutions are

for all t ≥ tfix ,

where Ffix is a reference to F(t) and tfix is a defined time. This limitation helps to determine steady-state values for the number of fertilized females at the end of the control process, which will promote stabilization in the number of individuals in the populations in the aquatic phase, as well as immature females and natural males. All nondominated solutions with fitness equal to 1 discovered by MGA_DENGUE are stored at each iteration (generation), forming an elite group. At the end of the optimization process, this elite will be the best approximation of the Pareto front found. In all iterations, part of the population (% pc) is copied to create an intermediate population (mating pool). The selection of individuals to be copied can be made through the Roulette Method. Roulette Method also known as Roulette Wheel Method, is used for selecting potentially useful solutions. Usually a proportion of the wheel is assigned to each of the solutions based on their fitness value. Then a random selection is made similar to how the roulette wheel is rotated and solutions with greater fitness values have better chance to be selected (more details can be found in Ref. [7]). In current algorithm greatest chance is given to individuals with fitness equal to 1. The individuals in the mating pool are used in the crossover operation. The genetic crossover operator allows individuals in the population to reproduce and form new generations. In MGA_DENGUE, crossover is effected through the random choice of the two individuals in an intermediate population and two crossover points (r1 and r2 ) (two point crossover) and, thus, the features of two solutions are recombined to form two offspring (solution 3 and solution 4). Among the four individuals (solution 1, solution 2, solution 3 and solution 4), two are chosen that have the best fitness for a return to the population, and these individuals replace the two worst. In the illustration below, we consider two individuals chosen for a crossover. Suppose two crossover points, r1 = 2 and r2 = 3, are obtained in a random process. Thus the two first elements in line 1 and the three first elements in line 2 of the solution 1 and solution 2 are separated. The elements of solution 1 and solution 2 recombine, generating solution 3 and solution 4, as shown in Fig. 5. After the crossover operator, individuals from the current generation are randomly chosen for mutation. The mutation probability should be set low (pm < 0.1). If the change in the value of the element is favorable, that is, if the number drawn is less than pm , then another drawing is performed to determine how many and which elements in each line of the solution shall undergo the change, and thereafter a new drawing takes place to determine the value of each of these elements.

r1=2

r2=3 Fig. 5. Illustration of the crossover operator on individuals.

82

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

Table 1 Parameters used for the implementation of the optimality problem (3). Parameter (day−1 ) ϕ C γ 0.5

13

0.07

r 0.5

β 1

βS

μA

0.7

0.05

μI 0.05

μF 0.05

Table 2 Parameters of the MGA_DENGUE.

μM 0.1

Parameters of the MGA_DENGUE T G n h

μS 0.1

120

The process is repeated G times and an approximate Pareto front is determined. The Pareto-optimal solutions are the solutions in elite. 5. Results We applied the MGA_DENGUE described to solve the optimization problem (3) with the following two objectives:

Minimize J1 =

T 0

Maximize J2 =

T

0

2

u1 + u22 + F 2 dt

(S2 )dt.

The objective here is to determine the u1 = u1 (t) and u2 = u2 (t) that minimize the function J1 , related to investments in insecticide, sterile mosquitoes and with the amount of fertilized females in the system (social cost) and to maximize J2 related to the preservation of the number of sterile mosquitoes introduced. The MGA_DENGUE was computationally implemented using MATLAB 7.6.0.324 (R2008a) [15] in a Dual Core i5-650 microcomputer with 4GB memory and 400GB hard drive in the Scientific Laboratory of Informatics (LCI, Laboratório Científico de Informática) in the Department of Biostatistics of the Institute of Biosciences at UNESP (Universidade Estadual Paulista) in Botucatu, São Paulo, Brazil (the average CPU time required is 10 min). Table 1 shows the parameters used for the implementation of the optimality problem. All parameters except r are measured in day−1 (per day). With the data in Table 1, we can calculate the initial conditions presented in (2) and the R value (R = 2.7778). Fig. 6 shows a simulation of the dynamic mosquito population without control for 300 days, starting with the following mosquito densities: A0 = 0.0001, I0 = 0, F0 = 0 and M0 = 0. In this situation, where there is a low density of eggs, larvae and pupae, the growth of

1 2 3 4

10

Density Density Density Density

of aquatic phase A(t) of immature female I(t) of fertilized females F(t) of natural male M(t)

300

10

K

pc

Pm

tfix

Ffix

50

0.70

0.05

60

0.30

where: T is the period of application of the control; G is the generation number; n is the number of individuals in the population; P is the penalty imposed upon unfeasible individuals; K is the number of elements in the elite set; pc is the probability of an individual of being selected for crossover; pm is the mutation probability; Ffix is the is a reference output on F(t); tfix is the defined time such that F(t) ࣘ Ffix for all t ࣙ tfix .

the population becomes accelerated from the 100th day. Stabilization starts from the 200th day. Table 2 presents the parameters used in the computational implementation of the MGA_DENGUE. The initial population was constructed using a random heuristic, where it created individuals from the population by random selections with u1 (t) ࢠ [0, 1] and u2 (t) ࢠ [0, 30] for all t ࢠ [0, T], with total period T equal to 120 days. The value T = 120 was chosen because it represents the time (in days) in which there is a greater infestation of dengue fever – hot and rainy periods in Brazil (about 4 months). The calculation of the integral was performed using Simpson’s 1/3rd Rule on discrete time varying from 0 to 120 days with unit spacing [12]. To solve the system of nonlinear ordinary equations (1) the Runge–Kutta 4th order method was used [12]. Using these methods and the parameters in Tables 1 and 2, the MGA_DENGUE was implemented to solve the optimization problem and the results are discussed as follows. Fig. 7 shows the Pareto front for the optimality problem. Each point in the Pareto front is associated with a curve from u1 and u2 . Taking three particular cases, as shown Fig. 7, we can show the combination of u1 and u2 that was chosen by MGA-DENGUE and the evolution of the mosquito population that these controls promoted (Fig. 8).

4.5 12

500

3

x 10

6

4

1

3.5

8

2.5

3

J2

Densities of mosquitoes

3

Case 3

2 6

1.5 4

1

4

Case 2

0.5 2

0 0

2 0

0

50

100

150

200

250

300

350

t

Fig. 6. The evolution of the mosquito population without control for 300 days, beginning with a low density of eggs, larvae and pupae.

Case 1 0.5

1

1.5

2

2.5 J1

3

3.5

4

4.5

5

x 10

4

Fig. 7. Illustration of Pareto front for the optimality problem associated with curves u1 and u2 in and three cases (solutions marked with “∗”). Case 1: small values in objectives inserted in J1 and J2 . Case 2: mid-sized values in objectives inserted in J1 and J2 . Large values in objectives inserted in J1 and J2 . The values in objectives J1 and J2 are mid-sized or large compared with Case 1.

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

0.35

30

Case 1 Case 2 Case 3

0.3

83

25

Case 1 Case 2 Case 3

3

0.2

Control variable u2

Control variable u1

0.25

1

0.15

0.1

2 0 0

2

15

2 10

1

1

0.05

20

5

3

20

40

60 Time (d)

80

100

2120

0 0

20

40

(a)

6 5 4 3

3

2

1

2 1

20

40

60 Time (d)

80

100

120

Evolution of the population immature females I

Evolution of the population in Aquatic phase A

7

0 0

120

Case 1 Case 2 Case 3

0.3

0.25

0.2

2

0.15

3 1

0.1

0.05

0 0

20

40

60 Time (d)

(c)

80

100

120

(d) 6 Case 1 Case 2 Case 3

2.5

2

1.5

2 1

3 1 20

40

60 Time (d)

80

100

120

(e)

Evolution of the population natural males M

3 Evolution of the population fertilized females F

100

0.35 Case 1 Case 2 Case 3

8

0 0

80

(b)

9

0.5

60 Time (d)

Case 1 Case 2 Case 3

5

4

3

2

2

3 1

0 0

1 20

40

60 Time (d)

80

100

120

(f)

Evolution of the population Machos Estereis S

300 Case 1 Case 2 Case 3

250

3

200

2

150

100

1

50

0 0

20

40

60 Time (d)

80

100

120

(g) Fig. 8. Solutions showed in Fig. 7 and defined as Case 1 (black), Case 2 (blue) and Case 3 (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

84

H.O. Florentino et al. / Mathematical Biosciences 258 (2014) 77–84

The functional J1 involves u1 (t), u2 (t) and F(t) which are variables related to amount of insecticide applied at time t, density of sterile male mosquitoes inserted in nature at time t and density female mosquitoes present in local at time t, respectively. Minimize the functional J1 implies in 3 aspects: determine the minimum amount of insecticide that should be used at each time t, i.e., minimize environmental damage that insecticide can cause on environment; determine the minimum amount sterile male mosquitoes to be inserted in nature (it could be an important decision in the process of production of these mosquitoes); u1 and u2 should be determined such that female mosquito population, F(t), decrease with the time t, because the female mosquitoes are directly responsible for the transmission of Dengue virus. The functional J2 is related to the amount of sterile male mosquitoes which remain in nature (preservation of sterile male mosquitoes) and is not desired that the insecticide acts on sterile male mosquitoes and reduce this population. Thus, a decision maker can choose solutions on the Pareto curve according to the interest in reduce J1 or increase J2 . Fig. 7 shows all efficient solutions proposed by MGA_DENGUE and we highlight three cases when small, middle and large values were obtained for objectives inserted into J1 and J2 , solutions marked with “∗”. We note in Fig. 8 that with combinations of u1 and u2 proposed by MGA_DENGUE, all segments of the population of natural mosquitoes decrease and the sterilized male population continues in greater amounts. The detached three cases indicate that must apply a larger amount of insecticide during 10 or 20 days, decreasing this amount after this period and then must insert sterile male mosquitoes until 70th day. In the next period should decrease the amount of sterile mosquitoes and insert little amount of insecticide, as shown in Fig. 8(a) and (b). Small values of u1 (t) means to apply little amount of insecticide at time t and big value of u2 (t) means to insert large amount of sterile male mosquitoes density at time t. Furthermore, a comparison of studied three cases suggests that when to need decrease the amount of sterile mosquitoes, should increase the amount of insecticide and vice versa, so that the natural population of mosquitoes continues decreasing over time. Another fact that Fig. 8 (from c to f) indicates that the control with insecticide is most effective in terms of lower the number of individuals in the population of mosquitoes than control with sterile males. But, on the other hand, the combination of these two controls has the advantage of reducing the use of insecticide and consequently reduces damage that it may cause in the environment, animals and human and do not induce insecticide resistance in the next generation of mosquitoes. In general, inserting constraints into optimal control problems complicates the resolution of differential equation systems. However, MGA-DENGUE offers great flexibility in the insertion of constraints, allowing a large gain in terms of control of the stated variables due to the insertion of the final conditions for F. A decrease can be observed in all segments of the mosquito population. The proposed multiobjective is a tool that offers a diversity of solutions, giving the decision maker an explanation of all possible optimized combinations using insecticide and sterile male that can be made. 6. Conclusion In the present work, we proposed numerical techniques using genetic algorithms to solve the model for multiobjective optimal control problems applied to the fight against the dengue mosquito.

In multiobjective, efficient solutions are calculated assuming different levels of investment in sterile male (J2 ), ranging in importance in order from small to large investments, giving decision makers the freedom to choose the most appropriate option for their investment interests, either in terms of economic costs or the social and environmental impacts of each control. In the present work, multiobjective optimization problems are exemplified with two performance functions, but MGA-DENGUE can be easily implemented with three or more functions, according to the user’s interest and without any change in its structure. The genetic algorithm is easy to implement and has been proved a good tool for determining optimum controls for the problem at hand. Acknowledgments We wish to thank FAPESP (Grant Nos. 2009/15098-0[Temático], 2010/07586-6 and 2014/01604-0), FUNDUNESP/PROPE/UNESP (Grant No. 0351/019/13), CNPq (Grant No. 303267/2011-9), CAPES and PROPG UNESP for their financial support. References [1] L.S. Barsante, R.T.N. Cardoso, J.L. Acebal, Otimização multiobjetivo no controle de gastos com inseticidas e machos estéreis no combate da dengue, in: Anais do XLIII Simpósio Brasileiro de Pesquisa Operacional, Ubatuba, SP, Brasil, 2011 Available at: http://www.xliiisbpo.iltc.br/pdf/87927.pdf (accessed 02.06.14). [2] L.S. Barsante, R.T.N. Cardoso, J.L. Acebal, Otimizando custos no combate a dengue através de algoritmo genético, in: X SBAI Simpósio Brasileiro de Automação Inteligente, São João Del Rei, MG Brasil, 2011 Available at: http://fei.edu.br/sbai/SBAI2011/87344.pdf (accessed 02.06.14). [3] A.C. Bartlett, Insect sterility, insect genetics and insect control, in: D. Pimentel (Ed.), Handbook of Pest Management in Agriculture, vol. II, CRC Press, Boca Raton, FL, 1990, pp. 279–287. [4] N. Becker, D. Petrick, C. Boase, J. Lane, M. Zgomba, C. Dahl, A. Kaiser, Mosquitoes and Their Control, 2nd ed., Springer, 2010. [5] M.A.L. Caetano, T. Yoneyama, Optimal and sub-optimal control in dengue epidemics, Optimal Control Appl. Methods 22 (2001) 63–73. [6] CDC – Centers for Disease Control and Prevention, 28 October 2010. Available at: http://www.cdc.gov/dengue (accessed 18.07.13). [7] K. Deb, Multi-objective Optimization using Evolutionary Algorithms, Wiley, Chichester, UK, 2001. [8] K. Dietz, The effect of immigration on genetic control, Theor. Popul. Biol. 9 (1976) 58–67. [9] L. Esteva, H.M. Yang, Mathematical model to assess the control of Aedes aegypti mosquitoes by the sterile insect technique, Math. Biosci. 198 (2005) 132–147. [10] L. Esteva, H.M. Yang, Control of dengue vector by sterile insect technique considering logistic recruitment, TEMA Tend. Mat. Apl. Comput. 7 (2006) 259–268. [11] C.M. Fonseca, P.J. Fleming, Multiobjective genetic algorithms, in: IEE Colloquium on ‘Genetic Algorithms for Control Systems Engineering’ (Digest No. 1993/130), 28 May 1993, IEE, London, UK, 1993. [12] C. Gerald, P. Wheatley, Applied Numerical Analysis, 7th ed., Addison-Wesley Publishing Company, 2003. [13] J.H. Holland, Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, MI, 1975. [14] E.F. Knipling, Sterile insect technique as a screwworm control measure: the concept and its development, in: O.H. Graham (Ed.), Symposium on Eradication of the Screwworm from the United States and Mexico, Misc. Publ. Entomol. Soc. America 62, College Park, MD, 1985, pp. 4–7. [15] MATLAB version 7.6.0.324 (R2008a), High Performance Numeric Computation and Visualization Software: Reference Guide, The MathWorks Inc., Natick, Massachusetts, 2008. [16] H. Pates, C. Curtis, Mosquito behavior and vector control, Annu. Rev. Entomol. 50 (2004) 53–70. [17] R.E. Plant, M. Mangel, Modeling and simulation in agricultural pest management, SIAM Rev. 29 (1987) 235–261. [18] H.S. Rodrigues, M.T.T. Monteiro, D.F.M. Torres, Dynamics of dengue epidemics when using optimal control, Math. Comput. Model. 52 (2010) 1667–1673. [19] J.D. Schaffer, Multiple objective optimization with vector evaluated genetic algorithms, in: Proceedings of the International Conference on Genetic Algorithm and their Applications 1985. [20] R.C.A. Thome, H.M. Yang, L. Esteva, Optimal control of Aedes aegypti mosquitoes by the sterile insect technique and insecticide, Math. Biosci. 223 (2010) 12–23. [21] WHO – World Health Organization. Dengue. Available at: http://www.who.int/ csr/resources/publications/dengue/001-11.pdf (accessed 19.07.13).