Original article Biomath 1 (2012), 1209262, 1–7 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus C. Dufourd ∗, Y. Dumont† ∗ IRD, CRVOI, Réunion Island, France Email: claire.dufourd@gmail.com † CIRAD, Umr AMAP, Montpellier, France Email: yves.dumont@cirad.fr Received: 10 July 2012, accepted: 26 September 2012, published: 31 December 2012 Abstract—To prevent epidemics of mosquito-transmitted diseases like Chikungunya in Réunion Island, we develop tools to control its principal vector, Aedes albopictus. Biological control tools, like the Sterile Insect Technique (SIT), are of great interest as an alternative to chemical control tools which are very detrimental to environment. The success of SIT is based on a good knowledge of the biology of the insect, but also on an accurate modeling of insects distribution. We model the mosquito dispersal with a system of coupled parabolic PDEs. Considering vector control scenarii, we show that the environment can have a strong influence on mosquito distribution and in the efficiency of vector control tools. Keywords-parabolic equation; existence; mosquito dis- persal; modeling; numerical simulation; splitting algo- rithm; vector control; Sterile Insect Technique. I. INTRODUCTION Chikungunya is an unusual vector-borne disease. First isolated in Tanzania in 1953, it is now geographically distributed in Africa, India and South-East Asia. After a huge epidemic in Réunion Island and in India in 2006, it appeared for the first time in Europe, in Italy, in 2007 (see [1], [2] and references therein). The symptoms that characterize Chikungunya are high fever, headache, persistant joint pain that can last several weeks. One way to reduce the risk of infection for the population is to control the vector populations. The principal vector for Chikungunya, is Aedes albopictus mosquito, commonly called the “Asian tiger” [3]. Standard vector control tools, like adulticide and lavicide, together with mechanical control are useful but cannot always be used for several reasons. Firstly, because Réunion Island is a hot spot of endemicity. Secondly, because mosquito can develop resistance to insecticides. Thirdly, because only approximately 10% of the island can be treated, due to its chaotic landscape. Therefore, it is necessary to consider new sustainable alternatives or additional tools, like the Sterile Insect Technique (SIT). SIT consists in releasing sterilized male mosquitoes that will mate with wild females which won’t be able to have offspring. Consequently, this will lead to the decrease of the vector population [4], [5]. The success of SIT is based on a good knowledge of the biology and the behavior of the vector, but also on an accurate modeling of its dispersal to optimize the impact of sterile males. The previous published models were temporal models, and didn’t take into account the spatial component. But, it is necessary to consider a spatio- temporal model to obtain realistic simulations and to simulate several vector control strategies. In a previous paper [6], we have considered a dispersal model with only adult females splitted in two compart- ments: the blood meal searchers and the breeding site searchers. This led to a system of two partial differential equations. Here, we add the aquatic stage, immature females, resting females, wild males, and finally, sterile males. This leads to a system of coupled (nonlinear) partial differential equations. After some theoretical re- sults and the presentation of the compartmental model, Citation: C. Dufourd, Y. Dumont, Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus, Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 1 of 7 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus we present briefly the numerical methods based on the operator splitting technique [7]. Finally, we end the paper with some numerical simulations with and without chemical or biological control. II. THE MATHEMATICAL MODEL Aedes albopictus is found in South-East Asia, the Pacific and the Indian Ocean islands, and up North through China and Japan. In the last twenty years, it has invaded several developed countries in Europe, the USA, Africa and South America (see [8] and references therein). When mosquitoes are not subject to stimuli, it is possible to assume that they move randomly in any direction [9]. This leads to a diffusion equation which can be extended to take into account the landscape heterogeneity or correlated random walk. Therefore, we have to incorporate advection terms or drift terms when mosquitoes, stimulated by attractants, move preferably in certain directions. For simplicity, we present the generic PDE that models the spread of a mosquito (sub)population. Of course, it is now well recognized that the environment heterogeneity can have an important effect [10], [11]. This is taken into account in the model by assuming spatial and temporal variations in the parameters. So, let u represent the number of insects. A possible modeling is to consider the following general advection-diffusion-reaction-like equation: ∂u ∂t = ∇ · (D∇u) − ∇ · (Cu) − ~V · ∇u − g (u) , (1) in a bounded domain Ω ⊂ Rn (where 1 ≤ n ≤ 3) with a smooth boundary ∂Ω. D (x, t) ≥ 0 is the diffusion (dispersion) coefficient or the diffusivity (D can even- tually depend on u). Entomologists usually admit that there is no passive transportation of Aedes albopictus mosquito by the wind. Conversely, the blood-seeking mosquitoes will follow odors and carbon dioxide carried by the wind; this is modeled by the term ~V (x, t) .∇u. Indeed, it is well known that carbon dioxide (CO2), in interaction with other components, acts as an attractant and induces a direct response to guide the mosquito towards the host. The breeding sites or the blood feeding sites attractions are modeled by the term ∇·(C (x, t) u). Thus it is necessary to take this important fact into account: the term ∇. · (C(x, t)u) represents a localized attractants (or a repellings) due to the presence or not of (blood or sugar) meals (like animals, humans or fruits...). C(x, t) represents an advection velocity toward favorable “places”. The definition of C takes into account wind’s direction and strength. For the sake of simplicity, the effective attraction areas are represented as ellipses. The attractor is set as one of the foci of the ellipse. The other focus point is calculated as a function of wind’s direction and strength. Outside the ellipse, there is no attraction and the related advection term is equal to zero. Note that if there is no wind, the attraction area is reduced to a disk of which the center is the attractor. The reaction term g can be nonlinear, and represents the time-evolution of the population (death-birth rates, migration, vector control...), and thus may depend on the mosquito population u, and some environmental parameters (temperature, position in the domain,...). We suppose that u (x, 0) = u0 (x) for x ∈ Ω, where u0 can be a continuous (or possibly discontinuous) function. It may be possible to consider generalized boundary conditions, like Robin conditions, −→ ∇u·~n + αu = ρ(x, t), for all (x, t) ∈ ∂Ω × (0, T ), where ~n stands for the exterior unit normal to the boundary ∂Ω. For the sake of simplicity we will consider homogeneous Neumann boundary conditions, i.e. α = 0. Finally, we deduce the following (quasi)linear parabolic equation:   ∂u (·, t) ∂t = ∇ · (D (·, t) ∇u) − ∇ (C (·, t) u) + ~V (·, t) · ∇u +g (u, ·, t) , x ∈ Ω and t > 0, u (x, y, 0) = u0 (x, y) ≥ 0, x ∈ Ω−→ ∇u · ~n = 0, for all x ∈ ∂Ω, and t > 0, (2) Problem (2) is a (quasi)linear parabolic equation, for which it is possible to show the existence of a local solution. Moreover, under mild conditions it may be possible to show that the solution is global [12]. A. The Compartmental Model In order to obtain some biologically interesting simu- lations, we take into account some biological facts about Aedes albopictus [3]. There are two main stages in the development of mosquitoes: an aquatic stage and an adult stage. The aquatic stage gathers eggs, larvae and pupae. The adult stage can be divided into several compart- ments: immature females, blood feeding females, breed- ing females, resting females and males. Since we assume no sex differences in the aquatic stage, mosquitoes, after emergence, are distributed between the immature female compartment and the male compartment, according to r, the ratio of adult females mosquitoes to the total mosquito population. After mating with males, immature females enter the feeding female compartment and seek Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 2 of 7 http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus ∂uA ∂t = NEgg(1 − uA K ) · 1b · ub(x, t)dt − (ηA + MA)uA, ∂uY ∂t = ∇ · (D∇uY ) + ~V · ∇uY − (MY + βY )uY + rηAuA, ∂uf ∂t = ∇ · (D∇uf ) + ~V · ∇uf − ∇(Cf (x)uf ) − (Mf + µf r1f )uf + µbf 1bub + βY uM uM + uM s uY when uM s 6= 0, = ∇ · (D∇uf ) + ~V · ∇uf − ∇(Cf (x)uf ) − (Mf + µf r1f )uf + µbf 1b · ub when uM s=0, ∂ur ∂t = ∇ · (D∇ur) + ~V · ∇ur − (Mf + µrb1b)ur + µf r · uf ,   ∂ub ∂t = ∇.(D∇ub) + ~V · ∇ub − ∇(Cb(x)ub) − (Mf + µbf 1b)ub + µrb · ur, ∂uM ∂t = ∇ · (D∇uM ) + ~V · ∇uM − ∇(Cf (x)uM ) − ∇(Cb(x)uM ) − Mm · uM + (1 − r)ηA · uA, ∂uM s ∂t = ∇ · (D∇uM s) + ~V · ∇uM s − ∇(Cf (x)uM s) − ∇(Cb(x)uM s) − Mms · uM s + Λs · 1s, ∀(x, t) ∈ QT uA(x, y, 0) = uA0 (x, y) x ∈ Ω uX (x, y, 0) = 0 x ∈ Ω with X ∈ {Y, f, b, M, M s} ~∇uX (x, y, t) · ~n = 0 x ∈ δΩ and 0 < t < T with X ∈ {A, Y, f, r, b, M, M s} (3) for blood meals before going into the resting compart- ment. Afterwards, the females pass into the breeding compartment seeking for a breeding site to deposit eggs. Once egg deposit is done, breeding females need blood and pass into the feeding compartment again. The eggs laid by the breeding females supply the aquatic stage. In order to simulate SIT control, we add a sterile male compartment. Sterile males are released in specific places. The transmission rate between the immature females and the feeding females is conditioned by the proportion of non-sterile males to the whole male popu- lation, near the immature females. Contrary to Anopheles mosquito, Aedes males are looking for young females and thus, in general, they are located near the breeding sites or near the hosts. We assume that resting females are not subjected to the attraction of blood meals or breeding sites. Resting females diffuse slowly, and their direction can be af- fected by the wind. The rate of transmission from one compartment to another allows us to take into account the average time spent by mosquitoes in each compart- ment. According to the previous explanations, we derive model (3). After rescaling, we consider Ω = [−1, 1]2, QT = Ω × (0, T ]. We set three indicator functions, 1b, 1f and 1s. Function 1b defines the area where the breeding females found a breeding site to lay eggs and become feeding females. 1f defines the area where feeding females found a blood meal and become resting females. 1s defines the area where sterile males are released. Cf represents the attraction due to blood meals, like houses, and Cb represents the attraction due to the breeding sites. MA, MY , Mf , Mm and Mms are respectively the mortality rates for the aquatic stage, the immature females, the mature females, the wild males and the sterilized males. ηA is the egg hatching rate, βY is the rate at which immature females become blood- feeding females, µf r is the rate at which blood-feeding females become resting females, µrb is the rate at which resting females become breeding females, µbf is the rate at which breeding females become blood-feeding females, NEgg is the average number of eggs laid per female, K is the carrying capacity of a breeding site, and Λs is the number of sterile males released periodically each τ days. Following [4], [5], we assume that the probability that an immature female becomes a feeding female depends on the ratio uM uM +uM s , when sterile males are released. System (3) can be rewritten in the following way, with u = (uA, uY , uf , ur, ub, uM , uM s) T ,   ∂u ∂t = ∇ · (D∇u) − ∇ (C (x) u) + ~V (x, t) · ∇u+ +Mu + γ(x, t) ∈ QT , ∇u · ~n = 0, x ∈ ∂Ω and 0 < t < T, u (x, 0) = u0 (x) x ∈ Ω, (4) Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 3 of 7 http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus with γ(x, t) = 0 BBBBBBB@ 0 0 0 0 0 0 Λs 1 CCCCCCCA , D (·, t) = D (·, t) 0 BBBBBBB@ 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 CCCCCCCA , C (·, t) = 0 BBBBBBB@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Cf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Cb 0 0 0 0 0 0 0 Cf + Cb 0 0 0 0 0 0 0 Cf + Cb 1 CCCCCCCA , M (·, t) = 0 BBBBBBB@ −(MA + ηA) 0 0 −rηA −(MY + βY ) 0 0 βY ( uM uM +uM s ) −(Mf + µf r1f ) 0 0 µf r1f 0 0 0 (1 − r)ηA 0 0 0 0 0 0 NEgg(1 − uA K )1b 0 0 0 0 0 0 0 0 0 0 −(Mf + µrb) 0 0 0 µrb −(Mf + µbf 1b) 0 0 0 0 −Mm 0 0 0 0 −Mms 1 CCCCCCCA . Problem (4) is a (quasi)linear weakly coupled par- abolic system, for which it is possible to show the existence of a local solution. The existence of a unique global solution may be proved using, for instance, [13]. To provide numerical simulations of system (3), we need to construct a reliable algorithm, that preserves most of the properties of the system (positivity of the solution, equilibrium, if any, and its (un)stability prop- erties, ...). B. The Numerical Methods We consider an operator splitting method [7]. This is an interesting method that enables us to solve separately each term of equation (1). So, we will solve successively the convective term, the diffusive term, and the reaction term, using the most efficient numerical method for each process. The full system can also be rewritten as follows ut = f (x, u) = A (u) + D (u) + R (u) , (5) where A represents the advective (or convective) terms, D the diffusive terms and R the reaction terms. Briefly, • the advection process is solved using the Corner Transport Upwind scheme (CTU) [14]. It is a total variation diminishing scheme, that preserves the monotonicity of the solution providing that we verify a CFL condition between the convective para- meters, the space-steps and the time-step. Moreover, this scheme minimizes the numerical diffusion and is even exact when ∆t is chosen appropriately. • the diffusion process is solved using the method of lines (MOL) for which we consider the second- order finite difference method for the spatial dis- cretization, and the TR-BDF2 (Trapezoide Rule - 2nd order Backward Difference Formula) method for the time discretization. • the reaction process is solved using the nonstandard finite difference method [15]. More details are developed in a forthcoming paper (see also [6]). Our scheme permits to provide several simulations with and without constant parameters. In particular, we show that the solution converges, at least numerically, to a steady state. The numerical algorithm is implemented in Scilab [16], while the visualization is obtained with “R” [17]. III. SIMULATIONS AND DISCUSSIONS We present some simulations in a homogeneous land- scape, with time independent parameters, where there are 4 breeding sites and 5 blood-feeding sites as in Fig. 1. We assume that all these attractors have the same attraction force and domains of attraction. For each attractor, the area of attraction is defined as an ellipse depending on the wind’s direction and strength. Note that, if there is no wind, the area of attraction is reduced to a disk. Fig. 1. Homogeneous landscape with 5 houses (blood meals) and 4 breading sites When the landscape is homogeneous, without wind, the mature females, i.e., the blood feeding females, the Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 4 of 7 http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus Fig. 2. Mature Females distribution with no wind, no vector control. Fig. 3. Mature Females distribution with North wind, mechanical control and releases of 1000 sterile males per week, from the 20th simulation day. resting females and the breeding females, tend to gather around the houses and the breeding sites (Fig. 2). This distribution is quite disturbed when wind and vector control get involved. Fig. 3, show the distribution of females when we consider North Wind, the removal of the two East breeding sites and the releases of sterile males. We notice that the majority of the mature females tend to migrate North, following odors or CO2 carried by the wind. They gather around the Western attractors since the two East breeding sites no longer exist. Moreover, we notice that the number of mature females has decreased. This decrease can be explained by the fact that sterile males have been released, but also because of the North migration that allows mosquitoes to leave the domain. SIT is respectful towards the environment, and is likely to be used as an alternative to the use of chemical products. On a temporal scale, we compare the impact of a 7-days periodic massive spraying of Deltamethrin around houses over a one-month treatment with 7-days Fig. 4. Mature Females distribution with North wind, and weekly releases of 1000 sterile males, near the North breading sites, from the 20th simulation day. Fig. 5. Mature Females distribution with North wind, and weekly releases of 1000 sterile males, near the South breading sites, from the 20th simulation day. pulsed sterile males releases near the breeding sites over the same period of time. Fig. 6 shows that with an appropriate choice in the periodicity of the releases, and the number of released sterile males, SIT could be a promising alternative to massive spraying. However, in order to have the most efficient results using SIT, it is important to have a good knowledge on the environment and take it into account. For instance, compare Fig. 4 and Fig. 5: they show the distribution of mature females controlled with SIT near the North breed- ing sites, and near the South breeding sites, respectively. We know that, with North wind, females tend to migrate towards the North; this is also the case for the released sterile males. That is why Northern releases have a lower impact on the number of mature females compared to Southern releases that are more efficient. Thus, it is more efficient to release sterile males downwind rather then Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 5 of 7 http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus Fig. 6. Comparison of the evolution of the number of mature females with a weekly massive spraying control, and SIT with weekly releases of 2000, 4000, 10000 and 20000 sterile males over one month. Fig. 7. Evolution of the density of mature females with Northern wind with two spatial strategies of SIT ( 1000 sterile males released per week) from the 20th day. upwind (on condition that the wind strength is not to high, otherwise mosquitoes will hide wherever they can rather than fly upwind). This result is even more obvious when we have a look at the temporal evolution of the number of mature females depending on which spatial strategy is chosen for SIT releases (Fig. 7). IV. CONCLUSION This work presents promising tool for modeling mos- quito dispersal. Thanks to the numerical simulations we could point out interesting results that can be a great use for optimizing SIT. Our work can be helpful to propose new field experiments; in particular it may be possible to consider temperature-varying parameters. This will be presented in a forthcoming paper. We could also add other compartments (sugar feeding compartment, sterile female compartment,... ). As we could see, the environment has a non-negligible influence on mosquito dispersal, and some works need to be done about land- scape ecology because, so far, little is known about the interactions between landscape, vegetation and Aedes albopictus dispersal [11]. ACKNOWLEDGMENT The TIS Project was financially supported by the French Ministry of Health and the FEDER Convergence Réunion 2007–2012 program. This paper is dedicated to Axel, for his fifteenth birthday. REFERENCES [1] Y. Dumont, F. Chiroleu, and C. Domerg, “On a temporal model for the chikungunya disease: Modeling, theory and numerics,” Mathematical Biosciences, vol. 213, no. 1, pp. 80–91 (2008). http://dx.doi.org/10.1016/j.mbs.2008.02.008 [2] Y. Dumont and F. Chiroleu, “Vector control for the chikungunya disease.” Mathematical Biosciences and Engineering, vol. 7, no. 2, pp. 313–345 (2010). http://dx.doi.org/10.3934/mbe.2010.7.313 [3] H. Delatte, C. Paupy, J. Dehecq, J. Thiria, A. Failloux, D. Fonte- nille et al., “Aedes albopictus, vector of chikungunya and dengue viruses in reunion island: biology and control,” Parasite (Paris, France), vol. 15, no. 1, pp. 3-13 (2008). [4] R. Anguelov, Y. Dumont, and J. Lubuma, “Mathematical mod- eling of sterile insect technology for control of anopheles mos- quito,” Computers and Mathematics with Applications, vol. 64, no. 3, pp. 374–389 (2012). http://dx.doi.org/10.1016/j.camwa.2012.02.068 [5] Y. Dumont and J. Tchuenche, “Mathematical studies on the sterile insect technique for the chikungunya disease and Aedes albopictus,” Journal of Mathematical Biology, vol. 65, no. 5, pp. 809–854 (2012). http://dx.doi.org/10.1007/s00285-011-0477-6 [6] Y. Dumont and C. Dufourd, “Spatio-temporal modeling of mosquito distribution,” in Proceedings of the 3rd Interna- tional Conference–AMiTaNS’11. AIP Conference Proceedings- American Institute of Physics, vol. 1404, pp. 162–165 (2011). [7] W. Hundsdorfer and J. Verwer, Numerical solution of time- dependent advection-diffusion-reaction equations. Springer Verlag, , vol. 33 (2007). [8] M. Benedict, R. Levine, W. Hawley, and L. Lounibos, “Spread of the tiger: global risk of invasion by the mosquito aedes albopictus,” Vector-Borne and Zoonotic Diseases, vol. 7, no. 1, pp. 76–85 (2007). http://dx.doi.org/10.1089/vbz.2006.0562 [9] P. Daykin, F. Kellogg, and R. Wright, “Host-finding and repul- sion of aedes aegypti,” The Canadian Entomologist, vol. 97, no. 3, pp. 239–263 (1965). http://dx.doi.org/10.4039/Ent97239-3 [10] G. Lemperiere, S. Boyer, and Y. Dumont, “Influence of rural landscape structures on the dispersal of the asian tiger mosquito Aedes albopictus : a study case at la reunion island,” 8th IALE (International Association of Landscape Ecology) World Congress 18–23 August 2011, Beijing, China. Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 6 of 7 http://dx.doi.org/10.1016/j.mbs.2008.02.008 http://dx.doi.org/10.3934/mbe.2010.7.313 http://dx.doi.org/10.1016/j.camwa.2012.02.068 http://dx.doi.org/10.1007/s00285-011-0477-6 http://dx.doi.org/10.1089/vbz.2006.0562 http://dx.doi.org/10.4039/Ent97239-3 http://dx.doi.org/10.11145/j.biomath.2012.09.262 C. Dufourd et al., Modeling and Simulations of Mosquito Dispersal. The Case of Aedes albopictus [11] Y. Dumont, “Modeling mosquito distribution. impact of the vegetation,” in Proceedings of ICNAAM 2011. AIP Conference Proceedings-American Institute of Physics, vol. 1389, pp. 1244– 1247 (2011). [12] A. Constantin and J. Escher, “Global solutions for quasilinear parabolic problems,” Journal of Evolution Equations, vol. 2, no. 1, pp. 97–111 (2002). http://dx.doi.org/10.1007/s00028-002-8081-2 [13] A. Constantin, J. Escher, and Z. Yin, “Global solutions for qua- silinear parabolic systems,” Journal of Differential Equations, vol. 197, no. 1, pp. 73–84 (2004). http://dx.doi.org/10.1016/S0022-0396(03)00165-7 [14] P. Colella, “Multidimensional upwind methods for hyperbolic conservation laws,” Journal of Computational Physics, vol. 87, no. 1, pp. 171–200 (1990). http://dx.doi.org/10.1016/0021-9991(90)90233-Q [15] R. Anguelov, Y. Dumont, and J. Lubuma, “On Nonstandard Finite Difference Method in Biosciences”, AMITANS 2012, AIP Conf. Proc. 1487, pp. 212–223 (2012). [16] Scilab Enterprises, Scilab: Le logiciel open source gratuit de calcul numrique, Scilab Enterprises, Orsay, France, 2012. [Online]. Available: http://www.scilab.org [17] R Development Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Com- puting, Vienna, Austria, 2011, ISBN 3-900051-07-0. [Online]. Available: http://www.R-project.org/ Biomath 1 (2012), 1209262, http://dx.doi.org/10.11145/j.biomath.2012.09.262 Page 7 of 7 http://dx.doi.org/10.1007/s00028-002-8081-2 http://dx.doi.org/10.1016/S0022-0396(03)00165-7 http://dx.doi.org/10.1016/0021-9991(90)90233-Q http://www.scilab.org http://www.R-project.org/ http://dx.doi.org/10.11145/j.biomath.2012.09.262 Introduction The Mathematical Model The Compartmental Model The Numerical Methods Simulations and Discussions Conclusion References