www.biomathforum.org/biomath/index.php/biomath REVIEW ARTICLE Optimal control strategies for a class of vector borne diseases, exemplified by a toy model for malaria Sebastian Aniţa∗, Edoardo Beretta†, Vincenzo Capasso‡ ∗Faculty of Mathematics, “Alexandru Ioan Cuza” University of Iaşi and “Octav Mayer” Institute of Mathematics of the Romanian Academy Iaşi, Romania sanita@uaic.ro †CIMAB, Interuniversity Centre for Mathematics Applied to Biology, Medicine, and Environmental Sciences, Italy edoardo.beretta@uniurb.it ‡ADAMSS, Centre for Advanced Applied Mathematical and Statistical Sciences Università degli Studi di Milano La Statale Milano, Italy vincenzo.capasso@unimi.it Received: 16 February 2019, accepted: 15 September 2019, published: 13 October 2019 Abstract—This paper contains a unified review of a set of previous papers by the same authors concerning the mathematical modelling and control of malaria epidemics. The presentation moves from a conceptual mathematical model of malaria trans- mission in an homogeneous population. Among the key epidemiological features of this model, two-age- classes (child and adult) and asymptomatic carriers have been included. As possible control measures, the extra mortality of mosquitoes due to the use of long-lasting treated mosquito nets (LLINs) and Indoor Residual Spraying (IRS) have been included. By taking advantage of the natural double time scale of the parasite and the human populations, it has been possible to provide interesting threshold results. In particular, key parameters have been identified such that below a threshold level, built on these parameters, the epidemic tends to extinction, while above another threshold level it tends to a nontrivial endemic state. The above model has motivated further analysis when a spatial struc- ture of the relevant populations is added. Inspired by the above, additional model reductions have been introduced, which make the resulting reaction- diffusion system mathematically affordable. Only the dynamics of the infected mosquitoes and of the infected humans has been included, so that a two-component reaction-diffusion system is finally Copyright: c©2019 Aniţa et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Sebastian Aniţa, Edoardo Beretta, Vincenzo Capasso, Optimal control strategies for a class of vector borne diseases, exemplified by a toy model for malaria, Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 1 of 19 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... taken. The spread of the disease is controlled by three actions (controls) implemented in a subdo- main of the habitat: killing mosquitoes, treating the infected humans and reducing the contact rate mosquitoes-humans. To start with, the problem of the eradicability of the disease is considered, while the cost of the controls is ignored. We prove that it is possible to decrease exponentially both the human and the vector infective population everywhere in the relevant habitat by acting only in a suitable subdomain. Later the regional control problem of reducing the total cost of the damages produced by the disease, of the controls and of the intervention in a certain subdomain is treated for the finite time horizon case. In order to take the logistic structure of the habitat into account the level set method is used as a key ingredient for describing the subregion of intervention. Here this subregion has been better characterized by both area and perimeter. The authors wish to stress that the target of this paper mainly is to attract the attention of the public health authorities towards an effective and affordable practice of implementation of possible control strategies. Keywords-malaria; nonlinear ODE models; qual- itative analysis; reaction-diffusion system; zero- stabilization; regional control; optimal control; epi- demic systems; I. INTRODUCTION ”Mathematical Models are lies that make us reach the truth”. (Paraphrased from Picasso) Human malaria is caused by one or a combina- tion of four species of plasmodia: Plasmodium falciparum, P. vivax, P. malariae, and P. ovale. The parasites are transmitted through the bite of infected female mosquitoes of the genus Anophe- les. Mosquitoes can become infected by feeding on the blood of infected people, and the parasites then undergo another phase of reproduction in the infected mosquitoes. According to [https://data.unicef.org/topic/ child-health/malaria/], “Prevention has been mainly carried out by the use of bed nets. Sleeping under insecticide-treated mosquito nets (ITNs) on a regular basis is one of the most effective ways to prevent malaria transmission and reduce malaria related deaths. Since 2000, production, procurement and delivery of ITNs, particularly Long Lasting Insecticide Treated Nets (LLINs) have accelerated, resulting in increased household ownership and use. ... Household ownership of ITNs/LLINs is uneven across countries in the region, with an average coverage of 66 per cent in sub-Saharan Africa ranging from less than 30 per cent to approximately 90 per cent but most countries have made considerable progress in the past decade. The proportion of sub-Saharan African households with at least one ITN increased to 67 per cent in 2016, thus, a third of households where ITNs are the main method of vector control did not have access to a net. Additionally, only 42 per cent of households had sufficient ITNs for all household members which is drastically short of the universal access of 100 per cent to this preventive measure.” This shows the relevance of studies as the ones proposed in this paper, concerning the identifica- tion of regions where universal access might be implemented, with respect to areas where this is not affordable because of various restrictions due to financial and logistic affordability. The earliest attempt to provide a quantitative un- derstanding of the dynamics of malaria transmis- sion was that of Ross [41]. Ross’ models consisted of a few differential equations to describe changes in densities of susceptible and infected people, and susceptible and infected mosquitoes. Macdonald [34], [35] extended Ross’ basic model, analyzed several factors contributing to malaria transmis- sion, and concluded that “the least influence is the size of the mosquito population, upon which the traditional attack has always been made”. The work of Macdonald had a very beneficial impact on the collection, analysis, and interpreta- tion of epidemic data on malaria infection [36] and guided the enormous global malaria-eradication campaign of his era. The classical Ross-Macdonald model has in- spired many contributions considering a variety of additional epidemiological features of the conta- gion [12], [26] (see also [13], [20], [22], [23], [27], Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 2 of 19 https://data.unicef.org/topic/child-health/malaria/ https://data.unicef.org/topic/child-health/malaria/ http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... [31], [42], [44]). A realistic mathematical model for an epidemic malaria system would require, at least for the human population, splitting the population with respect to stages of the disease, e. g. suscepti- bles, infectives, immunes, and alike; age structure, spatial structure, and alike. Even though we may suppose that indeed such a model is realistic, it might exhibit a complexity that make it hard for a sound mathematical analysis, and its validation as far as fitting real data is concerned. Expected uncertainty in the parameters would add additional complexity, so that the design of sound control strategies appears to be unaffordable. Reduction of complexity has then become a key issue, though one must take into account that reduction has not to lead to meaningless models with respect to real epidemic systems (a great interpretation of this concept has been offered by Picasso http://www.dailyartmagazine. com/pablo-picassos-bulls-road-simplicity/). This paper is meant to present a unified review of a set of previous papers by the same authors, the crucial issue being the implementation of optimal control strategies solely on an optimally chosen subregion of the habitat. In order to include specific issues of malaria epi- demic, we have started by considering a spatially homogeneous population with two age groups (children and adults) for the human population, and two classes of human individuals with respect to symptoms (asymptomatic and symptomatic) [39], [15]. In [39] the role of uncertainty in data was analyzed with respect to optimal control strategies. But the mathematical analysis of the full model was far from satisfactory, due to its intrinsic nontrivial structure. In [15], by taking advantage of the natural double time scale of the parasite and the human populations, the authors have been able to obtain a significant model reduction and then to provide interesting threshold results. In particular, there it is shown that key parameters can be identified such that below a threshold level the epidemic tends to extinction, while above another threshold level it tends to a nontrivial endemic state. In [8], [10] a spatial structure has been taken into account, but no age structure. Though the above mentioned simplifications lead de facto to a toy model, we have been using it as the basis for implementing sound realistic control strategies, inspired by the outcomes of the results in [39], [15]. The implementation of control strategies on an optimally chosen subregion of the habitat is a crucial addition of our proposal. Although we have tried to make clear the as- sumptions underlying our model, it has not been validated yet, by comparison with experimental data. Therefore it has to be cautioned that it is far from being the last word on a real malaria epidemic. However it is hoped that our model, with additional features that make it more realistic, when combined with efficient numerical methods, might provide the foundations for designing opti- mal control strategies by public health authorities. We wish to clarify that the reference to malaria does not exclude possible applications to other vector borne diseases, with relevant modifications that may take into account their specific epidemi- ological issues. For optimal control problems related to popu- lation dynamics we refer to the monographs [2], [3], [33], and references therein. A control problem of malaria for a space inde- pendent model may be found in [1]. For basic results and methods in shape optimiza- tion we refer to [25], [30], [40], [46]. It can be just mentioned that the models pre- sented here might be extended to other vector borne diseases, though this would require taking into account specific features of other diseases. This paper is organized as follows. In Section II a spatially homogeneous population with two age groups (children and adults) for the human population, and two classes of human individu- als with respect to symptoms (asymptomatic and symptomatic). By taking advantage of the double scale structure of the model, a reduced epidemic model is obtained, for which a qualitative analysis is carried out so to obtain a key threshold theorem. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 3 of 19 http://www.dailyartmagazine.com/pablo-picassos-bulls-road-simplicity/ http://www.dailyartmagazine.com/pablo-picassos-bulls-road-simplicity/ http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... Attention is paid to the behaviour of the total human population. In Section III the modelling, analysis and possible control strategies of the epidemic system with spatial structure has been reported. In particular a theorem is reported con- cerning the eventual eradicability of the disease, guaranteed by the proposed control strategies. In Section IV optimal regional control problems are discussed, including possible control strategies, and the optimal choice of the region of inter- vention. Different from previous papers, here the region of intervention has been characterized by both its area and perimeter. II. A COMPARTMENTAL MODEL WITH TWO AGE GROUPS The model developed in this paper combines two sub-models, the vectorial dynamics and the human host dynamics. Similar to Ross-MacDonald model [34], [32], once it is assumed that the vector population is comprised of only female Anopheles mosquitoes, they are categorized into Sv : uninfected (susceptible) mosquitoes, and Iv : infected mosquitoes. Therefore the total mosquito population is Nv = Sv + Iv. The total human population Nh has been di- vided into six classes: X : susceptible children , Y : susceptible adults, X1 : asymptomatic infected children, X2 :symptomatic infected children, Y1 : asymptomatic infected adults, Y2 : symptomatic infected adults (Y2). Hence, Nh = X + X1 + X2 + Y + Y1 + Y2. The parameters of the model are: ηh: Per capita birth rate of humans. ηv: Per capita birth rate of mosquitoes. λh: Transmission rate from infected mosquitoes to susceptible humans. λv: Transmission rate from infected humans to susceptible mosquitoes. a: Mosquito biting rate. b: Probability of transmission of infection from an infectious mosquito to a susceptible human. c: Probability of transmission of infection from an infectious human to a susceptible mosquitoes. fx: Proportion of symptomatic infection in children. fy: Proportion of symptomatic infection in adults. νx: Child natural recovery rate from asymptomatic infection. νy: Adult natural recovery rate from asymptomatic infection. rx: Child natural recovery rate from symptomatic infection. ry: Adult natural recovery rate from symptomatic infection. φ: Modification parameter for recovery rate from symptomatic infection. g: Human maturity rate. µv: Mosquito natural mortality rate. µh: Human natural mortality rate. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 4 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... δy: Disease induced mortality rate in adults. δx: Disease induced mortality rate in children. In the mathematical model we have denoted by λh the infection rate, due to the infective humans, acting on each susceptible mosquito, while by λv we have denoted the infection rate, due to the infective mosquitoes, acting on each susceptible human individual. We will adopt here the infection rates as in the classical Ross-McDonald model [34] so that λh = probability per unit time that a susceptible human is infected due a bite by an infective mosquito = relative abundance of mosquitoes with respect to the human population × probability per unit time a human is bitten by a mosquito × prob{the biting mosquito is infective| a mosquito has bitten a human} × prob{a susceptible human is infected| he has been bitten by an infective textnormalmosquito} = Nv Nh a Iv Nv b. Similarly, λv = probability per unit time that a susceptible mosquito is infected due to a bite to an infective human = probability per unit time a mosquito bites a human × prob{the bitten human is infective| a mosquito has bitten a human} × prob{a susceptible mosquito is infected| he has bitten an infective human} = a Ih Nh c. where Ih = X1 + X2 + Y1 + Y2. It is assumed that the individuals in the adult subpopulation have acquired an immunity and hence they are able to clear fast the infections (ry > rx), may have prolonged period of asymp- tomatic P. falciparum infection (1/νy > 1/νx) and suffers less disease induced mortality (δx > δy) compared to children. An extra death rate αv of the total mosquitoes population has been included, which may take into account the additional mortality due to the use of LLINs and the Indoor Residual Spraying (IRS) of mosquitoes. As in the force of infection of humans by mosquitoes, we have made the choice that this extra death rate of mosquitoes depends itself upon the relative abundance of mosquitoes with respect to the total human population. This modification takes into account the epidemiological issues of the Ross-McDonald model. In this framework it is more convenient to deal with fractional quantities by scaling each subpop- ulation with the total population at any time t as x(t) = X(t) Nh(t) , y(t) = Y (t) Nh(t) , · · · , iv(t) = Iv(t) Nv(t) . In this way, if we set m(t) := Nv(t) Nh(t) , the infection rates become λh(t) = bam(t)iv(t), λv(t) = acih(t). The compartmental mathematical model for malaria transmission is then represented by the system (1) of non-linear ordinary differential equa- tions, where x = 1 − (x1 + x2 + y + y1 + y2). Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 5 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... dx1 dt = (1 −fx)bamivx + (1 −φ)rxx2 − (g + νx)x1 −ηhx1 + x1(δxx2 + δyy2), dx2 dt = fxbamivx− (g + rx + δx)x2 −ηhx2 + δxx22 + δyx2y2, dy dt = gx + νyy1 + φryy2 − bamivy −ηhy + y(δxx2 + δyy2), dy1 dt = (1 −fy)bamivy + (1 −φ)ryy2 + gx1 −νyy1 −ηhy1 + y1(δxx2 + δyy2), (1) dy2 dt = fybamivy + gx2 − (ry + δy)y2 −ηhy2 + δxx2y2 + δyy22, div dt = λv − (λv + ηv)iv, sv = 1 − iv System (1) is complemented by the two equa- tions for the total human and mosquito populations dNh dt = (ηh −µh −δxx2 −δyy2)Nh, (2) dNv dt = (ηv −µv −αvm)Nv. (3) It is not difficult to show that System (1) is well posed, i.e., given suitable initial con- ditions, it admits a unique solution. More- over if the initial conditions are such that x(0),x1(0),x2(0),y(0),y1,y2(0) ≥ 0, with x(0) + x1(0) + x2(0) + y(0) + y1(0) + y2(0) = 1, and sv(0), iv(0) ≥ 0, with sv(0) +iv(0) = 1, then the same holds at any later time t > 0. We may then continue our analysis by refer- ring only to the lower dimensional vector xh := (x1,x2,y,y1,y2) T ∈ π := {xh ∈ R5+|x1 + x2 + y + y1 + y2 ∈ [0, 1]}, for the human population, and iv ∈ [0, 1] for the mosquito population. The total fraction of human susceptibles is sh = x + y. (4) If we introduce the vector of all fractions of infective humans ih := (x1,x2,y1,y2) T , (5) the total fraction of human infectives is given by ih = x1 + x2 + y1 + y2, (6) so that sh + ih = 1. (7) For the mosquito population we have just sv + iv = 1. (8) Our analysis may be concentrated on the time behaviour of the functions ih(t), and iv(t) de- scribing the fractions of total infective human population, and the fraction of infective mosquito population, respectively. A. The time scales According to the parameter values reported in [39], we can see that the natural birth rate ηh for the human population is of the order 10−5days−1, much lower to the natural birth rate ηv for the mosquito population which is of the order 10−2days−1, so that their ratio ε := ηh ηv � 1, say ε ∼ 10−3. It is then meaningful the introduction of two adimensional time scales as follows τh := ηht, τv := ηvt. The relationship between these two time scales is τh = ετv, so that τv, the natural time scale of the mosquito population results to be the fast time scale, while Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 6 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... τh, the natural time scale of the human population, is the slow time scale. By taking into account all the above, we can obtain the time evolution of the functions ih(t), and iv(t) with respect to these time scales. Let Ωh := {ih ∈ R 4 +| x1 + x2 + y1 + y2 ∈ [0, 1]} (9) and Ωv := [0, 1]. (10) Introduce the function f : Ωh × Ωv → R, (11) such that, for (ih, iv) ∈ Ωh × Ωv, f(ih, iv) := 1 ηh [F(ih, iv) −R(ih)], (12) with F(ih, iv) := bamiv(1 − ih) −ηhih, (13) and R(ih) := (1 − ih)(δxx2 + δyy2) +(νxx1 + νyy1) + φ(rxx2 + ryy2). (14) Because of (1), the time evolution of ih is then given by dih dτh = f(ih, iv), (15) coupled with div dτv = g(ih, iv), (16) where g(ih, iv) := 1 ηv [acih − (ηv + acih)iv]. (17) B. Qualitative analysis of the mathematical model At first we may provide some information about the time behaviour of m = Nv Nh . It can be shown that, at the fast time scale lim τv→+∞ Nv(τv) Nh(τv) = m∗, (18) for m∗ := ηv −µv αv = σv αv , (19) if we denote σv := ηv −µv > 0. (20) As a consequence, we may then assume, at the “slow” time scale τh, m(τh) := Nv(τh) Nh(τh) ≡ m∗. (21) 1) The slow time scale: With respect to the same slow time scale τh, equations (15) and (16) can be rewritten as follows dih dτh = f(ih, iv), ih(0) ∈ Ωh; (22) ε div dτh = g(ih, iv), iv(0) ∈ Ωv. (23) 2) The fast time scale: Viceversa, at the fast time scale the above equations become dih dτv = εf(ih, iv) (24) div dτv = g(ih, iv), (25) under the same initial conditions. Due to the fact that ε ' 0, at this scale ih can be considered constant, so that Equation (25) admits the equilibrium i∗v ≡ ϕ(ih) := acih ηv + acih , for any ih ∈ [0, 1]. (26) Notice that, for any ih ∈ [0, 1], ϕ(ih) is the only solution of g(ih, iv) = 0 Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 7 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... in [0, 1]. By the usual Lyapunov methods, the following proposition has been proven in [15]. Proposition II.1. For any choice of ih ∈ [0, 1], the point i∗v = ϕ(ih) is a globally asymptotically stable equilibrium solution for Equation (25). It is though interesting to directly evaluate the rate of convergence to the equilibrium i∗v = ϕ(ih) of Equation (25). If we take zv(τv) := iv(τv) − i∗v, (27) its evolution equation is dzv dτv = −K(ih)zv, (28) subject to zv(0) = iv(0) − i∗v. (29) Here K(ih) := ηv + acih ηv ≥ 1. (30) As a consequence, the solution of (28)-(29) is such that zv(τv) = zv(0) exp(−K(ih)τv) ≤ zv(0) exp(−τv), (31) i.e. |iv(τv) − i∗v| ≤ |iv(0) − i ∗ v|exp(−τv). (32) Equation (32) shows that an upper estimate of the convergence rate can be made, independent of ih. Furthermore, since |iv(0)−i∗v| ≤ 1, at the slow time scale τh this will be seen as |iv(τh) − i∗v| ≤ exp(−τh/ε). (33) which tends to zero extremely fast; indeed for ε ' 10−3, and τh = 10−1, exp(−τh/ε) ' exp(−100) ' 0. C. Qualitative analysis at the slow time scale Proposition II.1, and the fast convergence of iv(τh) to i∗v = ϕ(ih) at the slow time scale justify, by Tikhonov Theorem (see e.g.[14] ), the direct substitution of iv by i∗v in Equation (22). In addition, the above analysis allows us to assume that the ratio m(τh) := Nv(τh) Nh(τh) has reached its asymptotic value m∗ = σv αv , so that dih dτh = f̃(ih), ih(0) ∈ Ωh, (34) where f̃(ih) := 1 ηh [F̃(ih)) −R(ih)] (35) with F̃(ih) = −i2hα + ihβ ηv + acih . (36) Here α = ba2cm∗ + acηh, (37) and β = ba2cm∗ −ηhηv. (38) R(ih) is taken from (14). Due to the structure of (14), Equation (34) has always to be considered coupled with System (1) with (2) and (3), which makes its qualitative analysis rather complicated. A way to reduce such complexity is given below. We may observe that the function f̃ can be minorized by the function f`(ih) := 1 ηh [F̃(ih)) −Kih], (39) since R(ih) ≤ R(ih) := Kih, ih = ‖ih‖1 ∈ [0, 1], for K := max{δx + φry,νx} > 0. (40) On the other hand the function f̃ can be ma- jorized by the function fu(ih) := 1 ηh [F̃(ih)) −Kih], (41) since R(ih) ≥ R(ih) := Kih, ih = ‖ih‖1 ∈ [0, 1], Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 8 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... for K := min{νy,φrx} > 0. (42) We have denoted by ‖ih‖1 = x1 + x2 + y1 + y2. For the qualitative analysis it is then more convenient to study the minorant equation di`h dτh = f`(i`h), (43) together with the majorant equation diuh dτh = fu(iuh). (44) They are such that f`(ih) ≤ f(ih) ≤ fu(ih). (45) These equations are now decoupled from System (1). If we denote by ih(τh; ih(0)) the solution of Equation (34) ( coupled with System (1)-(3)), sub- ject to the inital condition ih(0); by i ` h(τh; ih(0)) the solution of (43), subject to the initial condi- tion ih(0) = ‖ih(0)‖1; and by iuh(τh; ih(0)) the solution of (44), subject to the initial condition ih(0) = ‖ih(0)‖1, then, by classical comparison theorems we may claim that for any τh ≥ 0, i`h(τh; ih(0)) ≤ ih(τh; ih(0)) ≤ i u h(τh; ih(0)). Then we have proved that, for any ih(0) ∈ (0, 1] : lim τh→+∞ [ i`h(τh; ih(0)), i u h(τh; ih(0)) ] = [̂ i`h, î u h ] where [̂ i`h, î u h ] is the ”attractor interval” for the total fraction ih(τh; ih(0)) of human infectives. D. The main threshold theorem By using the majorant equation, in [15] the authors have been able to show the following extinction theorem. Theorem II.2. If αv a2 ≥ bc ηv −µv ηv(ηh + K) (46) then the attractor interval [̂ ilh, î u h ] collapses to {0} and for all i.c. ih(0) ∈ (0, 1] we have: lim τh→+∞ ih(τh; ih(0)) = 0. (47) Since all parameters, but αv and a, are specific of the relevant populations, as expected we might control the malaria eradication acting on αv, the killing rate of mosquitoes by the use of LLIN’s and IRS, and on a, the mosquito biting rate, by the use of LLIN’s, and any other device/policy of contact reduction between mosquitoes and humans. We have to point out that Theorem II.2 offers only a sufficient condition for the eventual eradi- cation of the epidemic. By the minorant equation it has been also possi- ble to find a sufficient condition for the existence of a globally attractive endemic interval for the human infectives ih. E. About the total human population In the previous paragraphs we have reported a threshold theorem for the possible extinction of the malaria epidemic. An important issue is the persistence of the total human population, otherwise fractions may loose any epidemiological meaning. As a first step, let us rewrite Equation (2) at the slow time scale dNh dτh (τh) = σh − (δx x2(τh) + δy y2(τh)) ηh Nh(τh), (48) for τh ≥ 0, subject to an initial condition Nh(0) = Nh0 > 0. As above, we have denoted by σh := ηh−µh, that we assume to be strictly positive. In the sequel of this section we will write τ instead of τh. A minorant equation for (48) can be obtained as follows. If we denote by δ := max{δx,δy}, (49) and by ih(τ; ih0) the total fraction of human infectives, solution of Equation (34) (with (1)-(3)), Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 9 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... subject to an initial condition ih(0) ∈ Ωh, it is not difficult to realize that the following equation for a function N`(τ), τ ≥ 0, dN` dτ (τ) = δ ηh [ σh δ − ih(τ; ih0) ] N`(τ), (50) subject to the initial condition N`(0) = Nh0 > 0, is a minorant equation of (48) Hence, for any τ ≥ 0, the following inequality holds Nh(τ) ≥ N`(τ), (51) where N`(τ) :=Nh0 exp { δ ηh [ σh δ τ− ∫ τ 0 ih(u; ih0)du ]} . A sufficient condition for the “persistence in time” of the total human population Nh(τ), i.e. a sufficient condition which let us exclude that Nh(τ) ↓ 0+, as τ → +∞ (52) is (see [15]) σh > δ̂i u h, (53) where îuh is the upper extreme of the ”attractor interval”. Since in Theorem II.3 it occurs that î`h = î u h = 0 then we may claim that Theorem II.3. Under the assumptions of The- orem II.2 a sufficient condition for the persistence in time of Nh is σh > 0. (54) Unfortunately we have not been able to show that Nh is upper bounded, so that it might eventu- ally explode to +∞, though not in a finite time. Theorem II.2 suggests relevant control strate- gies for more detailed and realistic mathematical models. Next sections are devoted to the case of spatially structured models; eradication theorems will be shown, and an optimal control problem will be analyzed. III. A SPATIALLY STRUCTURED MODEL We now report about the modelling, analy- sis and possible control strategies of a malaria epidemic system with spatial structure, as from [8], [10]. About the modelling let us anticipate that, in order to keep the model affordable for the mathematical analysis further reductions have been introduced, still (we hope) in the spirit of Picasso. Indeed we have analyzed realistic control strategies, inspired by the above outcomes taken from [39], [15]. By considering a spatially structured system, we will refer to an habitat Ω ⊂ R2 (a nonempty bounded domain with a sufficiently smooth bound- ary ∂Ω). The two populations of humans and mosquitoes will be described in terms of their spatial densities. Specifically, we will consider a population of infected mosquitoes with spatial density u1(x,t) at a spatial location x ∈ Ω and a time t ≥ 0; while u2(x,t) will denote the spatial distribution of the human infective population. For our model reduction we have assumed that the spatial density C(x) of the total human popu- lation has been taken essentially constant in time, independent of removals by possible acquired im- munity, isolation, or death, so that C(x)−u2(x,t) provides the spatial distribution of susceptible hu- mans, at a spatial location x ∈ Ω and a time t ≥ 0. Further it has been assumed that the total susceptible mosquito population is so large that it can be considered time and space independent. As far as the “local incidence” for humans, at point x ∈ Ω, and time t ≥ 0, it is taken of the form (i.r.)H(x,t) = g(x,u1(x,t),u2(x,t)) = (C(x) −u2(x,t))h ( u1(x,t) C(x) ) , depending upon the local densities of both popula- tions via a suitable functional response h that will be better described later. Here we have taken into account that, in ac- cordance with the Ross-McDonald model [34] the function h, describing the force of infection of humans by mosquitoes, depends upon the relative Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 10 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... concentration of the total mosquito population with respect to the total human population, be- cause of the specific biting habits of humans by mosquitoes. On the other hand, we have here extended the (linear) response of the Ross-McDonald model, by using a possibly nonlinear functional response h( u1(x,t) C(x) ). This choice may allow possible satu- ration effects, σ−type responses, etc. (see [19], [17] ). For example behavioral changes can be taken into account; for a very large density of the infective population the force of infection represented by h(u1(x,t) C(x) ) may tend to reduce itself because of reduction of open exposure of the human population. This aspect had been already proposed in [19], which has motivated a large recent literature (see [28], and references therein). As pointed out in [36], seasonality of the ag- gressivity to humans by the mosquito population might also be considered in the functional response h; this has been a topic of the authors research programme in [18], [17], and [6]. As far as the “local incidence” for mosquitoes, at point x ∈ Ω, and time t ≥ 0, as in previous models [16], we assume that it is due to contagious bites to human infectives at any point x′ ∈ Ω of the habitat, within a spatial neighborhood of x represented by a suitable probability kernel k(x,x′), depending on the specific structure of the local ecosystem (see also [43]); as a trivial sim- plification one may assume k(x, ·) as a Gaussian density centered at x; hence the “local incidence” for mosquitoes, at point x ∈ Ω, and time t ≥ 0, is taken as (i.r.)M (x,t) = ∫ Ω k(x,x′)u2(x ′, t)dx′ . As proposed in [13], we have included spa- tial diffusion of the infective mosquito population (with constant diffusion coefficient to avoid purely technical complications), but we assume that the human population does not diffuse. More refined models may include some kind of mobility of the human population, due to migration mechanisms; but this would require further nontrivial complica- cies, that we leave to further investigations. Finally, we have denoted by η(x)u1(x,t) the (possibly space dependent) natural decay of the infected mosquito population, while a22u2(x,t) denotes the removal of the human infective popu- lation (by acquired immunity, isolation or death). All the above leads to the following over sim- plified model for the spatial spread of malaria epidemics, in which we have ignored various addi- tional features, such as the possible differentiation of the mosquito population, etc. (see e.g. [31])   ∂tu1(x,t) = d∆u1(x,t) −η(x)u1(x,t) + ∫ Ω k(x,x′)u2(x ′, t)dx′ ∂tu2(x,t) =−a22u2(x,t)+g(x,u1(x,t),u2(x,t)) in Ω × (0, +∞), where η(x) is the natural decay rate of the infected mosquitoes at x ∈ Ω (it is negative), a22 > 0, d > 0 are constants. The reader may recognize that the above system, apart from its spatial structure, corresponds to System (22)-(23), under additional simplifications. A. Regional Control Strategies The public health concern consists of providing methods for the eradication of the disease in the relevant population, as fast as possible. On the other hand, very often the entire domain Ω, of interest for the epidemic, is unknown, or diffi- cult to reach for the implementation of suitable environmental sanitation programmes. This has led the authors to suggest that implementation of such programmes might be done only in a given subregion ω ⊂ Ω, conveniently chosen so to lead to an effective eradication of the epidemic in the whole habitat Ω (“Think globally, act locally”, as in [4], [5], [6], [7]). This practice may have an enormous importance in real cases with respect to both financial and practical affordability. More recently, in [24], a patch control approach has been suggested, via the identification of zones with different disease burden. We assume here homogeneous Neumann bound- ary condition, to mean complete isolation of the Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 11 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... habitat: ∂νu1(x,t) = 0 on ∂Ω × (0, +∞), where ∂ν denotes the normal derivative. We also impose initial conditions{ u1(x, 0) = u 0 1(x) u2(x, 0) = u 0 2(x) in Ω, where u01 and u 0 2 are the initial densities of the population of infective mosquitoes and of the human infected population, respectively. Let ω ⊂ Ω be a nonempty open subset; we denote by χω the characteristic function of ω and use the convention χω(x)w(x) = 0, x ∈ R2 \ω, even if function w is not defined on the whole set R2 \ω. Our goal is to study the controlled system  ∂tu1(x,t) −d∆u1(x,t) = η(x)u1(x,t) + ∫ Ω k(x,x′)u2(x ′, t)dx′ −γ1(x,t)χω(x)u1(x,t), ∂tu2(x,t) =−(a22 +γ2(x,t)χω(x))u2(x,t) +(1−γ3(x,t)χω(x))g(x,u1(x,t),u2(x,t)), (55) for (x,t) ∈ Q, together with ∂νu1(x,t) = 0, (x,t) ∈ Σ (56) and u1(x, 0) = u 0 1(x), u2(x, 0) = u 0 2(x),x ∈ Ω, (57) where Q = Ω × (0, +∞), Σ = ∂Ω × (0, +∞). The idea underlying the controls is the follow- ing: γ1(x,t)u1(x,t)χω(x) represents the additional killing of mosquitoes by the combined action of a general insecticide spray- ing, and the use of treated nets in the subregion ω; γ2(x,t)u2(x,t)χω(x) describes the treatment of the infected population in the subregion ω; γ3(x,t)χω(x) is the reduction of the contact rate mosquitoes- humans by means of treated nets. As far as the epidemiological significance of the various control terms, let us consider at first the quantity γ1(x,t). This may be interpreted as a harvesting rate of infective mosquitoes, by the use of chemical-physical devices such as insecticides, and γ1(x,t)u1(x,t) as the corresponding distroyed mosquitoes population at location x ∈ Ω, and time t > 0. The control γ2(x,t) represents the recovery rate due to the medical treatment at location x ∈ Ω, and time t > 0; the corresponding cured population is γ2(x,t)u2(x,t). The control γ3(x,t) gives the additional seg- regation rate for the population of infective mosquitoes and for the human susceptible popu- lation at location x ∈ Ω, and time t > 0, due to the use of the treated bed nets (see e.g. [45], [47], and references therein). 1) Wellposedness: We have been working un- der the following technical assumptions: (H1) η ∈ L∞(Ω), k ∈ L∞(Ω × Ω), k(x,x′) ≥ 0 a.e. in Ω × Ω,∫ Ω k(x,x′)dx > 0 a.e. x′ ∈ Ω; (H2) C ∈ L∞(Ω), C(x) ≥ τ > 0 a.e. in Ω, where τ is a positive constant; (H3) h : R → [0, +∞), and h|[0,+∞) is continuously differentiable and increasing, h(x) = 0, for x ∈ (−∞, 0], h(x) ≤ a21x, for any x ∈ [0, +∞), where a21 is a positive constant; (H4) u01, u 0 2 ∈ L ∞(Ω), 0 ≤ u01(x), 0 ≤ u02(x) ≤ C(x) a.e. x ∈ Ω. Based on epidemiological considerations ([10]), the control functions γ = (γ1,γ2,γ3), introduced above, have been assumed to belong to G = {γ = (γ1,γ2,γ3)∈L∞(Q)×L∞loc(Q)×L ∞(Q) : 0 ≤ γ1(x,t) < Γ1, 0 ≤ γ2(x,t), 0 ≤ γ3(x,t) < Γ3, a.e. (x,t) ∈ Q}, for suitable constants Γ1, Γ3 > 0. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 12 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... Theorem III.1. [Existence and Uniqueness] For any γ ∈G, Problem (55) has a unique (weak) solution (uγ1,u γ 2 ). The solution satisfies 0≤uγ1 (x,t), 0≤u γ 2 (x,t)≤C(x), a.e.(x,t)∈Q. 2) Eradicability: Definition III.2. We say that the disease is eradicable if for any u10 and u 2 0 satisfying As- sumption (H4) there exists a γ ∈ G such that the solution (uγ1,u γ 2 ) to the controlled system satisfies lim t→+∞ u γ 1 (·, t) = lim t→+∞ u γ 2 (·, t) = 0 in L∞(Ω). In [10] it has been shown that Theorem III.3. For any u10,u 2 0 satisfying hy- pothesis (H4) and such that u10(x),u 2 0(x) > 0 a.e. x ∈ Ω \ω, and for any γ = (γ1,γ2,γ3) ∈ G, the solution (uγ1,u γ 2 ) satisfies u γ 1 (·, t) → 0 and u γ 2 (·, t) → 0 in L∞(Ω) (as t → +∞) slower than or as slow as exp (−a22t). The above theorem provides information on the lower bound for the convergence to zero of the epidemic. Additional analysis based on the magnitude of the principal eigenvalue of an ap- propriate operator, provides information about an upper bound on the convergence. For a given nonempty open subregion ω of Ω, such that Ω\ω 6= ∅, a given constant control γ = (γ1,γ2,γ3) ∈ G, and ε ∈ (0,a22), consider the following eigenvalue problem   −d∆ψ1(x) − (a22 −ε)ψ1(x) −η(x)ψ1(x) − ∫ Ω k(x,x′)ψ2(x ′)dx′ + γ1χω(x)ψ1(x) = λψ1(x),x ∈ Ω ∂νψ1(x) = 0,x ∈ ∂Ω (ε + γ2χω(x))ψ2(x) −(1 −γ3χω(x))a21ψ1(x) = 0,x ∈ Ω, By the Krein-Rutman Theorem, we may state that Proposition III.4. The principal eigenvalue λω1ε of the above problem is real, simple and admits at least an eigenvector (ψ1,ψ2) such that ψ1(x) ≥ ψ01, ψ2(x) ≥ ψ 0 2, a.e. x ∈ Ω, for suitable positive constants ψ01,ψ 0 2. By usual comparison theorems, the following monotonicity properties can be obtained. Proposition III.5. For a given set of control parameters γ ∈G, (i) The mapping ε ∈ (0,a22) 7→ λω1ε is continu- ous and strictly increasing. (ii) The mapping ω ⊂ Ω 7→ λω1ε is increasing, with respect to the ordering by inclusion. (iii) If {ωn}n∈N∗ is an increasing sequence of open subsets of Ω such that ∪∞n=1ωn = Ω, then lim n→+∞ λωn1ε = λ Ω 1ε. So that we may finally state our main eradica- bility result, i.e. an exponential eradicability. Theorem III.6. If ε ∈ (0,a22), γ ∈ G, and ω ⊂ Ω are such that the principal eigenvalue of the above problem is λω1ε > 0, then lim t→+∞ u γ 1 (·, t) = lim t→+∞ u γ 2 (·, t) = 0, in L∞(Ω) faster than or as fast as exp(−(a22 − ε)t). The following remark guarantees that this theo- rem is indeed applicable. Remark III.7. If Γ1 > a22 + ‖η‖L∞(Ω), then for ε ∈ (0,a22) sufficiently small, γ1 sufficiently close to Γ1, γ2 sufficiently large and ω sufficiently large, we have that λω1ε ≥ 0, and consequently the constant control γ = (γ1,γ2,γ3) (γ3 is arbitrary but fixed in [0, Γ3)) can make the epidemic expo- nentially eradicable. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 13 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... IV. OPTIMAL REGIONAL CONTROL OF THE EPIDEMIC For an optimal control problem during a finite time interval of intervention [0,T], we need to introduce a cost functional which takes into ac- count all relevant costs concerning on one side the costs of intervention associated with the functions (γ1,γ2,γ3) ∈ G; on the other side they must be taken also into account the costs deriving from loss of work hours, hospitalization, and alike asso- ciated with the infected human population (`(u2)), and the general costs associated with the specific choice of the subregion of intervention ω ⊂ Ω ( by assigning a suitable function α(x), x ∈ Ω, the specific costs related to the logistic structure of the habitat can be taken into account) . It should be then of the form J(γ,ω) = ∫ T 0 ∫ Ω ζ1(γ1(x,t))C(x)χω(x)dxdt + ∫ T 0 ∫ Ω ζ2(γ2(x,t)u2(x,t))χω(x)dxdt + ∫ T 0 ∫ Ω ζ3(γ3(x,t))C(x)χω(x)dxdt + ∫ T 0 ∫ Ω `(u2(x,t))dx dt + ∫ Ω α(x)χω(x)dx. (58) Based on suitable epidemiological considera- tions (see [10]), the following choices can be made for the specific cost functionals ζ1,ζ2,ζ3: (H5) ζ1 : [0, Γ1) → [0, +∞) is a continuously differentiable function, bijective, strictly in- creasing; ζ3 : [0, Γ3) → [0, +∞) is a continuously differentiable function, bijective, strictly in- creasing; (H6) ζ2 : [0, +∞) → (Γ2, Γ̃2] is a continuous differentiable function, bijective, and strictly decreasing. A. Optimal Control for a Fixed Region of Inter- vention At first let the subregion ω ⊂ Ω be fixed, so that the cost functional J(γ,ω) is just a function of γ, say J(γ). We shall denote by GT = {γ = (γ1,γ2,γ3)∈L∞(QT )3 : 0 ≤ γ1(x,t) < Γ1, 0 ≤ γ2(x,t), 0 ≤ γ3(x,t) < Γ3 a.e. (x,t) ∈ QT}, where QT = Ω × (0,T). Under the above assumptions the following the- orem holds. Theorem IV.1. The optimal control problem Minimizeγ∈GT J(γ) admits at least one solution γ∗ = (γ∗1,γ ∗ 2,γ ∗ 3 ). If we denote by dJ(γ) the gradient of J(γ), an optimal choice γ∗ is such that dJ(γ∗)(w) ≥ 0, for any w as in next proposition. Proposition IV.2. For any γ ∈ GT and w ∈ L∞(QT )×L∞(QT )×L∞(QT ) such that γ+θw ∈ GT for any θ > 0 sufficiently small, we have that dJ(γ)(w) = ∫ T 0 ∫ Ω w1(x,t)u γ 1(x,t)q1(x,t)χω(x)dxdt + ∫ T 0 ∫ Ω w1(x,t)C(x)ζ ′ 1(γ1(x,t))χω(x)dxdt + ∫ T 0 ∫ Ω w2(x,t)ζ ′ 2(γ2(x,t)u γ 2(x,t))u γ 2(x,t)χω(x)dxdt + ∫ T 0 ∫ Ω w2(x,t)u γ 2 (x,t)q2(x,t)χω(x)dxdt + ∫ T 0 ∫ Ω w3(x,t)C(x)ζ ′ 3(γ3(x,t))χω(x)dxdt + ∫ T 0 ∫ Ω (w3(x,t)(C(x) −u γ 2 (x,t)) ×h( u γ 1 (x,t) C(x) )q2(x,t)χω(x))dxdt, where (q1,q2) is the solution to (59) Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 14 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ...   ∂tq1(x,t) + d∆q1(x,t) = −η(x)q1(x,t) + γ1(x,t)χω(x)q1(x,t) −(1 −γ3(x,t)χω(x)) C(x) −uγ2 (x,t) C(x) h′ ( u γ 1 (x,t) C(x) ) q2(x,t), ∂tq2(x,t) = − ∫ Ω k(x′,x)q1(x ′, t)dx′ + (a22 + γ2(x,t)χω(x))q2(x,t) +(1−γ3(x,t)χω(x))h ( u γ 1 (x,t) C(x) ) q2(x,t)+ζ̃ ′ 2(γ2(x,t)u γ 2 (x,t))γ2(x,t)χω(x)+` ′(u γ 2 (x,t)), for (x,t) ∈ QT ; ∂νq1(x,t) = 0, (x,t) ∈ ΣT ; q1(x,T) = q2(x,T) = 0, x ∈ Ω. (59) Numerical results have been obtained by the gradient method, for suitable choices of all param- eters (see [10]). B. Optimal Choice of the Region of Intervention The optimization with respect to both the con- trol functions and the subregion of intervention by the gradient method, requires the evaluation of the directional derivative of the cost functional with respect to both the parameters γ1,γ2,γ3, and the region ω. A convenient way to handle the shape and position of ω is to use the level set method [40]. According to the implicit representation of subsets of Ω, we consider those subsets ω for which there exists a smooth function ϕ : Ω → R such that ω = {x ∈ Ω; ϕ(x) > 0} and ∂ω = {x ∈ Ω; ϕ(x) = 0}. Hence, instead of investigating the total cost function J defined in the first section, we shall deal with (H denotes the usual Heaviside function) J̃(γ,ϕ) = ∫ T 0 ∫ Ω ζ1(γ1(x,t))C(x)H(ϕ(x))dxdt + ∫ T 0 ∫ Ω ζ2(γ2(x,t)u2(x,t))H(ϕ(x))dxdt + ∫ T 0 ∫ Ω ζ3(γ3(x,t))C(x)H(ϕ(x))dxdt + ∫ T 0 ∫ Ω `(u2(x,t))dxdt + ∫ Ω α(x)H(ϕ(x))(x)dx. (60) where now γ ∈ GT , ϕ : Ω → R is continuously differentiable, and (u1,u2) is the solution to   ∂tu1(x,t) −d∆u1(x,t) = −η(x)u1(x,t) + ∫ Ω k(x,x′)u2(x ′, t)dx′ −γ1(x,t)H(ϕ(x))u1(x,t), ∂tu2(x,t) =−(a22 +γ2(x,t)H(ϕ(x)))u2(x,t) +(1−γ3(x,t)H(ϕ(x)))g(x,u1(x,t),u2(x,t)), (61) for (x,t) ∈ QT , subject to the boundary and initial conditions ∂νu1(x,t) = 0, (x,t) ∈ ΣT (62) u1(x, 0) =u 0 1(x), u2(x, 0) =u 0 2(x), x∈Ω, (63) where QT = Ω × (0,T), ΣT = ∂Ω × (0,T). Actually H is usually replaced by its mollified version Hε(s) = 12 (1 + 2 π arctan s ε ) (ε > 0 is a sufficiently small number). It is possible to prove that [10] Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 15 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... Theorem IV.3. For any γ ∈GT and w ∈ L∞(QT ) ×L∞(QT ) ×L∞(QT ) such that γ + θw ∈GT for any θ > 0 sufficiently small, and for any smooth functions ϕ,ψ : Ω → R we have that dJ̃ε(γ,ϕ)(w,ψ) = ∫ T 0 ∫ Ω w1(x,t)Hε(ϕ(x))u γ,ϕ 1 (x,t)q1(x,t)dxdt + ∫ T 0 ∫ Ω w1(x,t)C(x)ζ ′ 1(γ1(x,t))Hε(ϕ(x))dxdt + ∫ T 0 ∫ Ω w2(x,t)ζ ′ 2(γ2(x,t)u γ,ϕ 2 (x,t))u γ,ϕ 2 (x,t)Hε(ϕ(x))dxdt + ∫ T 0 ∫ Ω w2(x,t)u γ,ϕ 2 (x,t)q2(x,t)Hε(ϕ(x))dxdt + ∫ T 0 ∫ Ω w3(x,t)C(x)ζ ′ 3(γ3(x,t))Hε(ϕ(x))dxdt + ∫ T 0 ∫ Ω w3(x,t)(C(x) −u γ,ϕ 2 (x,t))Hε(ϕ(x))h ( u γ,ϕ 1 (x,t) C(x) ) q2(x,t)]dxdt + ∫ Ω ψ(x)[δε(ϕ(x))Ψ(x)]dx, (64) where δε(x) = H ′ ε(x) = 1 π ε x2 + ε2 , and Ψ(x) is given by Ψ(x) = ∫ T 0 [ C(x)ζ1(γ1(x,t))+ζ2(γ2(x,t)u γ,ϕ 2 (x,t))+C(x)ζ3(γ3(x,t))+α(x)+γ1(x,t)u γ,ϕ 1 (x,t)q1(x,t) +γ2(x,t)u γ,ϕ 2 (x,t)q2(x,t) + γ3(x,t)(C(x) −u γ,ϕ 2 (x,t))h ( u γ,ϕ 1 (x,t) C(x) ) q2(x,t) ] dt. Here (uγ,ϕ1 ,u γ,ϕ 2 ) is the solution to (61), and (q1,q2) is the solution to  ∂tq1(x,t) + d∆q1(x,t) = −η(x)q1(x,t) + γ1(x,t)Hε(ϕ(x))q1(x,t) −(1 −γ3(x,t)Hε(ϕ(x)) C(x) −uγ,ϕ2 (x,t) C(x) h′( u γ,ϕ 1 (x,t) C(x) )q2(x,t), ∂tq2(x,t) = − ∫ Ω k(x′,x)q1(x ′, t)dx′ + (a22 + γ2(x,t)Hε(ϕ(x))q2(x,t) +(1−γ3(x,t)Hε(ϕ(x))h ( u γ,ϕ 1 (x,t) C(x) ) q2(x,t)+ζ̃ ′ 2(γ2(x,t)u γ,ϕ 2 (x,t))γ2(x,t)Hε(ϕ(x)) + ` ′(u γ,ϕ 2 (x,t)), for (x,t) ∈ QT ; ∂νq1(x,t) = 0, (x,t) ∈ ΣT ; q1(x,T) = q2(x,T) = 0, x ∈ Ω. (65) Remark 1: We may notice that for a given γ = (γ1,γ2,γ3) ∈GT the quantity ∂sϕ(x,s) = −δε(ϕ(x,s))Ψ(x) (66) for x ∈ Ω, s > 0, where (u1,u2) is the solution of the state equation and (q1,q2) is the solution to the adjoint equation as above, in terms of the fictitious ”time” variable s, gives the descent of the optimal search with respect to ϕ. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 16 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... We may introduce additional costs due to shape of the subregion ω, by considering the perimeter of that region. As in [29] and [9], this would add the following term in the cost functional J̃(γ,ϕ) (see (60)) ∫ Ω β(x)δ(ϕ(x)) |∇ϕ(x)|dx. (67) As a consequence, in the gradient dJ̃ε(γ,ϕ)(w,ψ) (see (64)) last term would then become∫ Ω ψ(x)[δε(ϕ(x))Ψ̂(x)]dx + ∫ ∂Ω ψ(x) δε(ϕ(x)) |∇ϕ(x)| ∂νϕ(x)d` (68) where now Ψ̂(x) = ∫ T 0 [C(x)ζ1(γ1(x,t)) + ζ2(γ2(x,t)u γ,ϕ 2 (x,t)) +C(x)ζ3(γ3(x,t)) + γ1(x,t)u γ,ϕ 1 (x,t)q1(x,t) +γ2(x,t)u γ,ϕ 2 (x,t)q2(x,t) +γ3(x,t)(C(x)−u γ,ϕ 2 (x,t))h ( u γ,ϕ 1 (x,t) C(x) ) q2(x,t) +α(x) −β(x)div ( ∇ϕ(x) |∇ϕ(x)| ) ]dt. Under these circumstances Equation (66) becomes ∂sϕ(x,s) = −δε(ϕ(x,s))Ψ̂(x), x ∈ Ω,s > 0, (69) subject to the boundary condition β(x) δε(ϕ(x,s)) |∇ϕ(x)| ∂νϕ(x) = 0, x ∈ ∂Ω,s > 0 (70) (see [29] and [9]). We may recognize that this is strictly related to the literature on image segmentation (see [38], [37], [21]). V. CONCLUSIONS • Based on the main ideas of relevant literature concerning epidemiological issues of malaria, we have proposed a mathematical model de- scribing the dynamics of infected mosquitoes and humans in a spatially structured habitat. • In order to reduce the number of infected mosquitoes and humans, we have taken into account three possible control measures to be implemented only in a suitable subdomain ω of the relevant global habitat. • To start with, we have shown that if such a subdomain of intervention is sufficiently large and if the magnitude of the control efforts is sufficiently large, then eventual eradication of both infected populations is possible at an exponential rate. • We have then analyzed the optimal control problem for the reduction of the infected human population (in a finite time interval) at a minimum cost. • Further analysis for the identification of an optimal subdomain ω, together with optimal control efforts, has been left to future inves- tigations. — ACKNOWLEDGMENT This work has been carried out in the framework of the European COST project CA16227: “In- vestigation and Mathematical Analysis of Avant- garde Disease Control via Mosquito Nano-Tech- Repellents”. The authors would like to thank the anonymous referees for their valuable comments. REFERENCES [1] Agusto FB, Marcus N, Okosun KO. Application of op- timal control to the epidemiology of malaria. Electron. J. Diff. Eqs. 2012;81:1–22. [2] Aniţa S. Analysis and Control of Age-Dependent Popu- lation Dynamics. Dordrecht: Kluwer Acad. Publ.; 2000. [3] Aniţa S, Arnăutu V, Capasso V. An Introduction to Opti- mal Control Problems in Life Sciences and Economics. From Mathematical Models to Numerical Simulation with MATLAB. Basel: Birkhäuser; 2011. [4] Aniţa S, Capasso V. A stabilizability problem for a reaction-diffusion system modelling a class of spatially structured epidemic systems. Nonlin. Anal. Real World Appl. 2002;3:453–464. [5] Aniţa S, Capasso V. A stabilization strategy for a reaction-diffusion system modelling a class of spatially structured epidemic systems (think globally, act locally). Nonlin. Anal. Real World Appl. 2009;10:2026–2035. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 17 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... [6] Aniţa S, Capasso V, On the stabilization of reaction-diffusion systems modeling a class of man-environment epidemics: a review. Math. Meth. Appl. Sci. 2010;33:1235–1244. [7] Aniţa S, Capasso V. Stabilization of a reaction-diffusion system modelling a class of spatially structured epi- demic systems via feedback control. Nonlin. Anal. Real World Appl. 2012;13:725–735. [8] Aniţa S, Capasso V. Stabilization of a reaction-diffusion system modelling malaria transmission. Discrete Cont. Dyn. Syst. Ser. B, Special issue in honour of “A. Fried- man” 2012;17:1673–1684. [9] Aniţa S, Capasso V, Dimitriu G. Controlling an alien predator population by regional controls. Nonlin. Anal. Real World Appl. 2019; 46: 82–97. [10] Aniţa S, Capasso V, Dimitriu G. Regional Control for a spatially structured malaria model. Math. Meth. Appl. Sci. 2019;42:2909–2933. [11] Arnăutu V, Neittaanmäki P. Optimal control from theory to computer programs. Dordrecht: Kluwer Academic Publ.; 2003. [12] Aron JL, May RM. The population dynamics of malaria. In: “Population Dynamics of Infectious Diseases. The- ory and Applications” (ed. R.M. Anderson). London: Chapman & Hall 1982:139–179. [13] Bacaer N, Sokhna C. A reaction-diffusion system mod- eling the spread of resistance to an antimalarial drug. Math. Biosci. Eng. 2005;2:227–238. [14] Banasiak J, Lachowitz M Methods of Small Parameter in Mathematical Biology, Boston: Birkhauser; 2014. [15] Beretta, E, Capasso, V, Garao D.G. A mathematical model for malaria transmission with asymptomatic car- riers and two age groups in the human population. Mathematical Biosciences 2018;300:87–101. Erratum. Mathematical Biosciences 2018;303:155–156. [16] Capasso V. Asymptotic stability for an integro- differential reaction-diffusion system. J. Math. Anal. Appl. 1984;103:575–588. [17] Capasso V. Mathematical Structures of Epidemic Sys- tems, 2nd revised printing, Lecture Notes Biomath., Vol. 97. Heidelberg: Springer-Verlag; 2009. [18] Capasso V, Maddalena L. Periodic solutions for a reaction-diffusion system modelling the spread of a class of epidemics. SIAM J. Appl. Math. 1983;43:417–427. [19] Capasso V, Serio G. A generalization of the Kermack- McKendrick deterministic epidemic model. Math. Biosci. 1978;42:43–61. [20] Chamchod F, Britton NF. Analysis of a vector-bias model on malaria transmission. Bull. Math. Biol. 2011;73:639–657. [21] Chan TF, Vese LA. Active contours without edges IEEE Transactions on Image Processing 2001;10:266–277. [22] Chitnis N, Cushing JM, Hyman JM. Bifurcation analysis of a mathematical model for malaria transmission. SIAM J. Appl. Math. 2006;67:24–45. [23] Chitnis N, Smith TA, Steketee RW. A mathematical model for the dynamics of malaria in mosquitoes feed- ing on a heterogeneous host population. J. Biol. Dyn. 2008;2:259–285. [24] Chitnis N, Shapira A, Shindler C. Mathematical analysis to prioritise strategies for malaria elimination. J. Theor. Biol. 2018;455:118–130. [25] Delfour MC, Zolesio JP. Shapes and Geometries. Met- rics, Analysis, Differential Calculus and Optimization. 2nd Edition. SIAM, Philadelphia; 2011. [26] Dietz K. Mathematical models for transmission and control of malaria. In “Principles and Practice of Malar- iology” (eds. W. Wernsdorfer, Y. McGregor). Churchill Livingstone, Edinburgh; 1988:1091–1133. [27] Dietz K, Molineaux T, Thomas A. A malaria model tested in the African savannah. Bull. World Health Org. 1974;50:347–357. [28] d’Onofrio A, Manfredi P. Information-related changes in contact patterns may trigger oscillations in the en- demic prevalence of infectious diseases. J. Theor. Biol. 2009;256:473–478. [29] Getreuer P. Chan-Vese segmentation. Image Processing On Line 2012;2:214–224. [30] Henrot A, Pierre M. Variation et Optimisation de Formes, Mathématiques et Applications. Berlin: Springer-Verlag, Berlin; 2005. [31] Killeen GF, McKenzie FF, Foy BD, Schieffelin C, Billingsley PF, Beier JC. A simplified model for pre- dicting malaria entomological inoculation rates based on entomologic and parasitologic parameters relevant to control. Am. J. Trop. Med. Hyg. 2000;62:535–544. [32] Kingsolver, J G Mosquito host choice and the epidemi- ology of malaria. American Naturalist 1987;130:811– 827. [33] Lenhart S, Workman JT. Optimal Control Applied to Biological Models. Chapman and Hall; 2007. [34] McDonald G. The analysis of equilibria in malaria. Trop. Dis. Bull. 1952;49:813–829. [35] Macdonald G. The Epidemiology and Control of Malaria. London: Oxford University Press; 1957. [36] Molineaux L, Gramiccia G. The Garki Project. Research on the Epidemiology and Control of Malaria in the Sudan Savanna of West Africa. Geneva: World Health Org.; 1980. [37] Morel JM, Solimini S. Variational Models for Image Segmentation: with seven image experiments. Boston: Birkhauser; 1994. [38] Mumford D, Shah J. Optimal by piecewise smooth functions and associated variational problems. Comm. Pure Appl. Math. 1989;42:577–685. [39] Mwanga, G, Haario, H, Capasso, V. Optimal control problems of epidemic systems with parameter uncertain- ties: Application to a malaria two-age-classes transmis- sion model with asymptomatic carriers Mathematical Biosciences 2015;261:1–12. [40] Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. New York: Springer; 2003. [41] Ross R. The Prevention of Malaria. 2nd edition. Murray, London; 1911. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 18 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 S. Aniţa, E. Beretta, V. Capasso, Optimal control strategies for a class of vector borne diseases ... [42] Ruan S, Xiaob D, Beierc JC. On the delayed Ross- Macdonald model for malaria transmission. Bull. Math. Biol. 2008;70:1098–1114. [43] Shcherbacheva A, Haario H, Killeen GF. Modeling host- seeking behavior of African malaria vector mosquitoes in the presence of long-lasting insecticidal nets. Math. Biosci. 2018;295:36–47. [44] Smith T, Killeen G, Maire N, Ross A, Molineaux L, Tediosi F, Hutton G, Utzinger J, Dietz K, M. Tanner M. Mathematical modeling of the impact of malaria vac- cines on the clinical epidemiology and natural history of Plasmodium Falciparum malaria: overview. Am. J. Trop. Med. Hyg. 2006;75(2):1–10. [45] Sochantha T, Hewitt S, Nguon C, Okell L, Alexan- der N, Yeung S, Vannara H, Rowland M, Socheat D. Insecticide-treated bednets for the prevention of Plas- modium falciparum malaria in Cambodia: a cluster- randomized trial. Trop. Med. Int. Health 2006;11:1166– 1177. [46] Sokolowski J, Zolesio JP, Introduction to Shape Opti- mization. Berlin: Springer-Verlag; 1992. [47] White LJ, Maude RJ, Pongtavornpinyo W, Saralamba S, Aguas R, Van Effelterre T, Day NPJ, White NJ. The role of simple mathematical models in malaria elimination strategy design. Malaria J. 2009;8:212. [48] WHO Malaria Report; World Health Organization; 2017. Biomath 8 (2019), 1909157, http://dx.doi.org/10.11145/j.biomath.2019.09.157 Page 19 of 19 http://dx.doi.org/10.11145/j.biomath.2019.09.157 Introduction A compartmental model with two age groups The time scales Qualitative analysis of the mathematical model The slow time scale The fast time scale Qualitative analysis at the slow time scale The main threshold theorem About the total human population A Spatially Structured Model Regional Control Strategies Wellposedness Eradicability Optimal Regional Control of the Epidemic Optimal Control for a Fixed Region of Intervention Optimal Choice of the Region of Intervention CONCLUSIONS References