Regional Control for Spatially Structured Mosquito Borne Epidemics

This paper is a continuation of a previous paper by the first two authors to appear in the same IWR Special Issue for Scientific Computing. We are concerned with an optimal regional control problem for spatially structured vector borne epidemic system, considering malaria as a case study. A conceptual reduced mathematical model of malaria had been presented consisting of a two-component reaction-diffusion system. Three actions (controls) had been included: killing mosquitoes, treating the infected humans and reducing the contact rate mosquitoes-humans. The problem which is faced concerns the optimal choice of the region of intervention, by introducing a cost functional which takes into account the total cost of the damages produced by the disease, of the controls and of the intervention in a certain subdomain, for a finite time horizon case. A gradient algorithm had been proposed for the search of a minimal value of the cost functional, with respect to both the control parameters and the region of intervention. The scope of the present paper concerns the numerical implementation of such an algorithm. The level set method has played a major role for the mathematical description of the subregion of intervention. The outcomes of a series of numerical simulations are reported, under a variety of parameter scenarios.


Introduction
This paper contains a continuation of [2] in which an optimal control problem had been considered for a reaction diffusion system modelling vector borne epidemics, with malaria as a working example.
By considering a spatially structured system, we refer to an habitat ⊂ R 2 (a nonempty bounded domain with a sufficiently smooth boundary ∂ ).
The two populations of humans and mosquitoes are described in terms of their spatial densities. Specifically, u 1 (x, t) denotes the spatial density of the population of infected mosquitoes at a spatial location x ∈ and a time t ≥ 0; while u 2 (x, t) denotes the spatial distribution of the human infective population. The spatial density C(x) of the total human population has been assumed constant in time, so that C(x) − u 2 (x, t) provides the spatial distribution of susceptible humans, at a spatial location x ∈ and a time t ≥ 0. It has been assumed that the total susceptible mosquito population is so large that it can be considered time and space independent.
The proposal by the authors, going back to a series of papers [1,3], consists of reducing the area of intervention to an optimally chosen subregion ω of the entire habitat .
As anticipated in [3] and [2], the controlled system is the following one.
Let ω ⊂ be a nonempty open subset; we denote by χ ω the characteristic function of ω and use the convention even if the function w is not defined on the whole set R 2 \ ω.
The control γ 2 (t) represents the healing rate due to the medical treatment at time t > 0, so that γ 2 (t)u 2 (x, t)χ ω (x) describes the treatment of the infected population in the subregion ω.
Finally, the control γ 3 (t) gives the additional segregation rate for the population of infective mosquitoes and for the human susceptible population at time t > 0, due to the use of the treated bed nets, so that γ 3 (t)χ ω (x) is the reduction of the contact rate mosquitoes-humans by means of treated nets (see e.g. [11,12], and the references therein).
As a consequence of the above, we shall assume that Since we are going to face a finite time horizon for the control, we will refer to the following set for the control parameters As concerns the costs for the controls, in [3] we have proposed specific models in order to take into account realistic situations concerning specific epidemic cases. For malaria control specific cost functions ζ 1 , ζ 2 , ζ 3 have been assumed satisfying the following assumptions.
Here is the plan of the paper. In Section 2 the regional control problem is presented, paying attention to the mathematical description of the relevant subregion by the level set method. The gradient of the cost functional obtained in [2] with respect to both the controls and the subregion of intervention has been reported. This has driven the computational issues reported in Section 3 together with results from numerical simulations within a variety of parameter scenarios. damages produced by the disease and of the intervention in the subset ω in the time interval [0, T ] are considered, then we get the total cost function: Here γ ∈ G T and ω is a subset of . As far as the last three terms in (1) are concerned, the term T 0 l(u 2 (x, t))dx dt is meant to take into account the costs, over the whole domain and the total observation/control period, deriving from loss of work hours, hospitalization, drugs, etc.; while the terms α area (ω) + β perimeter (ω) take into account costs of transport of intervention devices (nets, chemicals, personnel, etc.) which may depend on the geometry of the subregion of intervention ω. For our scopes, the geometry of a planar set is indeed characterized by its area and perimeter; the coefficients α, and β take into account the specific logistic structure of the habitat; we may assume α, β ≥ 0.
We shall treat the problem of reducing the value of J with respect to γ and ω. From a computational point of view a convenient way to handle the shape and position of the subregion ω is the level set method (see [9], and the references therein). This is carried out by the implicit representation of ω. A convenient way to handle the shape and position of ω is to use the implicit representation according to 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 above, we shall deal withJ Here H is the Heaviside function. To describe the perimeter of the region ω, we have referred to the Modica-Mortola formula [7] (ε is a sufficiently small parameter). We shall approximate this cost function by the following one, where H is replaced by its mollified version H σ (s) = 1 2 1 + 2 π arctan s σ (σ > 0 is a sufficiently small number) Here u γ,ϕ (2) A rigorous proof of the existence of a global optimum of the above cost functional is a very hard task. A usual search algorithm is based on the gradient method. The following theorem provides the gradient with respect to both γ and φ; its proof can be found in [2].
for any θ > 0 sufficiently small, and for any smooth functions ϕ, ψ : → R we have that and Here q γ,ϕ Based on the above theorem, the descent towards an optimal γ can be obtained by considering and truncating if necessary.
The descent towards an optimal ϕ can be obtained by introducing an artificial time θ. We introduce a functionφ(x, θ) such that, given an initial guess ϕ old (x), we takeφ(x, 0) = ϕ old (x). At each iteration as ϕ new (x), we take the solution to the following PDE subject to the boundary condition and the initial conditionφ (x, 0) = ϕ old (x).
we may take, as boundary condition, instead of (4),

Computational Issues
In the first paragraph, we describe the numerical methods adopted for the solution of systems (2) and (3). In the second paragraph, we provide the values of all parameters used in the numerical tests. In the last paragraph, we report and discuss the computational results. The optimal control procedure implemented here is based on the gradient method presented in Section 2 and the consequent conceptual algorithm described in [2]. The reported numerical simulations have been performed with an in-house MATLAB code. We wish to stress here that the search of the optimal subregion ω is a very hard computational problem; mathematically the search should explore a very highly dimensional space of shapes. This is a well known issue in the realm of statistical shape analysis [5]. Such an intrinsic complexity has been reduced via the choice of the initial function φ within a parameterized family of shapes (e.g. ellipses), hence in a finite dimensional space; then the search goes on by looking for the parameters of the chosen shape that optimize the relevant cost functional. This may explain why in the reported numerical simulations the final optimal subset ω keeps the same shape as the one given via the initial function φ. From a purely mathematical point of view the search is supposed to move among all possible shape classes, which is clearly unaffordable from a computational point of view.

Numerical Methods
Systems (2) and (3) are discretized by finite elements in space and finite differences in time.

Space Discretization
We first apply a standard Galerkin procedure to the weak formulations of systems (2) and (3). The unit square domain is discretized by a uniform grid of 80 × 80 bi-linear finite elements (Q1), yielding a total amount of 6561 discretization nodes. The stiffness matrix is computed exactly, whereas the mass matrix is obtained by applying the mass lumping technique.

Time Discretization
After space discretization by finite elements, we obtain the semidiscrete problem that consists of a system of ordinary differential equations (ODEs). We solve these ODEs by employing a first order semi-implicit finite difference scheme, where the linear diffusion and reaction terms are approximated by Backward Euler, whereas the non-linear reaction terms are approximated by Forward Euler. As a result, at each time step it is required the solution of a linear system of algebraic equations of dimension 2 × 6561 degrees of freedom. The linear system is solved by the Gaussian elimination with the built-in function of Matlab. The time step size is t = 0.05.
For further details on the numerical approximation of parabolic equations we refer e.g. to [10].

Parameter Calibration
Except otherwise stated, the values of the parameters are given as detailed here below.
Functions η(x 1 , x 2 ) and C(x 1 , x 2 ), representing the growth rate of the infected mosquitoes and the spatial density of the total human population, respectively, are given by η(x 1 , x 2 ) = −10  and Functions ζ 1 , ζ 2 and ζ 3 are given by (see also Fig. 2): with α 1 = 0.1, 1 = 500 and b 1 = 1; with α 2 = 0.01, 2 = 500, a 2 = 0.02 and c 2 = 600; with α 3 = 0.01, 3 = 0.85 and b 3 = 1.  Functions h and l are taken as h(s) = 4s and l(s) = s, whereas constants α, β, d and a 22 are set to α = 100, β = 0.01, d = 1e − 4 and a 22 = 0.5. We recall that α is the coefficient in front of the cost functional term H (ϕ), which penalizes large values of the area of the control region. Regarding the diffusion coefficient d, its small value (1e − 4) yields a reaction-diffusion problem with dominant reaction. However, using higher diffusion Table 5 Test 2, optimal control, Case 4. Evolution of cost functional, control parameters γ 1 , γ 2 , γ 3 and number of infected mosquitos (# mosq) and humans (# hum) at final simulation time T = 1 as a function of the iterations (iter) of the optimal control procedure iter γ 1 γ 2 γ 3 cost # mosq # hum coefficients values, we believe that non-realistic diffusion effects might occur, with a too fast spread of the outbreak. The initial distributions of infected mosquitos and humans, used in all next numerical simulations, are plotted in Fig. 1 and are set as follows (Fig. 2):

Test 1: Reduction of the Infected Mosquitoes/Humans Distributions
In this first test, we keep fixed the control region ω with and we investigate the influence of the control parameters γ 1 , γ 2 and γ 3 on the cost functional and on the evolution of the infected mosquitos and humans distributions. The kernel function k(x, x ) is the Dirac function centered in x, thus The results reported in Table 1 and Fig. 3 show that the most sensitive parameter is γ 2 , whose increase yields a significant decay of the number of infected mosquitos and humans, with respect to non-controlled system.

Test 2: Optimal Control
In this test, we apply the optimal control procedure in four different scenarios (Cases 1-4), detailed here below, depending on the choice of the kernel function k(x, x ) and initial shape of function ϕ.
-Case 1: -Case4:ϕ(x 1 , x 2 ) = 0 has the quasi-elliptical shape reported in Fig. 7 (first row, second columns) and k x, x = δ x − x . For each of the four cases, we perform two optimal control procedures: one minimizing the cost functional only with respect to the parameters γ 1 , γ 2 and γ 3 (called in the following (γ 1 , γ 2 , γ 3 ) control strategy); and one minimizing the cost functional both with respect to the parameters γ 1 , γ 2 and γ 3 and to the function ϕ (called in the following (γ 1 , γ 2 , γ 3 , ϕ) control strategy).
The results reported in Tables 2, 3 -except Case 1, the minimal cost is achieved by performing the (γ 1 , γ 2 , γ 3 , ϕ) control strategy; -in all four cases, the area of the optimal control region is significantly reduced with respect to the initial control region.
We observe that the optimization process is highly sensitive with respect to the area of the control region, which is always significantly reduced in the optimal configuration with respect to the initial one, whereas the shape and location of the optimal region remain similar to the initial one. Thus, the simulations and the optimal control procedures show that acting only on a subregion, with a significantly smaller area than the entire domain, might control successfully the epidemics.
Moreover, the results show that, when the control region is located at the center of the domain (Cases 1, 2, 4), where the epidemics is stronger, the optimal configuration has a much smaller area than in Case 3, where the control region is shifted with respect to the center of the epidemics.

Future Directions
A difficult task which has emerged during the current research on malaria has been the identification of parameters for a real epidemic. For the case of malaria, some parameters have been identified in the literature, as reported in [4,8]. In a near future, we would like to implement the algorithms as from the above, with a full set of required parameters.

Concluding Remarks
As in our previous literature, first of all we have shown, supported by numerical simulations in a variety of parameter scenarios, that it is indeed possible to eventually eradicate a manenvironment epidemic, by acting only on the environment, while anyhow medically treating the human infected population. As stressed in [2], we do not claim that the simplified model presented here as a working example contains all details which might be required to act on a real epidemic system; still our hope is that our simulations may convince the public health authorities that it is possible, and it can be affordable, to control the full system by implementing the control strategies only in a suitably chosen subregion of the relevant habitat.
It is worth reporting here a statement (abridged) taken from [6] concerning the possible role of our investigations: "a model is only an approximate interpretation of reality and it is always wrong in some small or relevant elements. The destiny of the model presented here is to be rapidly improved thanks to novel knowledge coming from new observations and better assumptions. The Authors hope that many and more brilliant minds will read the present pages, will identify and highlight putative mistakes, will get inspiration for their research, and will produce better, more complete, and useful models. . . . . . .. If the speculations presented here on implications for surveillance, control, and therapy of [malaria] will contribute, even only minimally, to save . . . .. human life. . . . . . . . . , then the Authors have accomplished their small mission."