Optimal control for disease vector management in SIT models: an integrodifference equation approach

Vector-borne diseases are a major public health concern inflicting high levels of disease morbidity and mortality. Vector control is one of the principal methods available to manage infectious disease burden. One approach, releasing modified vectors (such as sterile or GM mosquitoes) Into the wild population has been suggested as an effective method of vector control. However, the effects of dispersal and the spatial distribution of disease vectors (such as mosquitoes) remain poorly studied. Here, we develop a novel mathematical framework using an integrodifference equation (discrete in time and continuous in space) approach to understand the impact of releasing sterile insects into the wild population in a spatially explicit environment. We prove that an optimal release strategy exists and show how it may be characterized by defining a sensitivity variable and an adjoint system. Using simulations, we show that the optimal strategy depends on the spatially varying carrying capacity of the environment.


Introduction
Vector borne diseases are a major public health concern, causing high levels of morbidity and leading to nearly one million deaths, annually (WHO 2016). Often vector control methods are the only feasible management option with the aim of vector control being to suppress or often eliminate insect vector populations. The sterile insect technique (SIT) is a well established empirical method for reducing population size. This technique has been widely used to suppress or eradicate many insect pest species, B Michael B. Bonsall michael.bonsall@zoo.ox.ac.uk 1 Mathematical Ecology Research Group, Department of Zoology, University of Oxford, Oxford OX1 3PS, UK through the release of modified insects (lab-sterilised, irradiated and/or novel genetic technologies) . Control is achieved as mating between sterile and wildtype mosquitoes reduces the reproductive potential of the wild population.
In evaluating the efficacy of SIT, mathematical modelling is an effective tool (Manoranjan and van den Driessche 1986;Isidoro et al. 2009;Yakob and Bonsall 2009;Legros et al. 2012;Li and Zou 2012;Potgieter et al. 2013;Li and Yuan 2015). SIT is an area-wide control method and hence involves spatial domains. Different continuous or discrete spatial model types have been proposed to investigate the effects of space on the efficacy of SIT. The choice between different model formulations depends on the specific biological problem (see Table 1). An early example is that developed by Manoranjan and van den Driessche (1986) who analysed the effect of non-uniform sterile insect releases in a heterogeneous environment. Another modelling approach explored considered random releases across the environment and concluded that the effectiveness of SIT is highly influenced by spatial heterogeneity (Ferreira et al. 2008). More recently, Yakob and Bonsall (2009) developed a spatially explicit model where the release of sterile insects was assumed to be uniform across space. Their research suggests that high levels of insect dispersal reduces the effectiveness of the SIT. Legros et al. (2012) analysed and compared two release strategies: (1) spatially uniform releases and (2) releases at discrete locations. They concluded that uniform releases are more effective than releasing at specific locations.
However, a constant dilemma is whether model predictions replicate real biological situations as the field data are often incomplete or not straightforward to analyse. In order to achieve successful vector control, we need a much more thorough understanding of how mosquitoes would respond to these control interventions (such as sterile insect releases) by focusing more on their biology, ecology and behaviour . Spatial spread is a key element in mosquito reproduction and some studies argue that SIT technique may not be very effective in controlling mosquito populations due to their dispersal and distribution (Ferreira et al. 2008).
To address this issue we present an alternative modelling framework that allows us to consider population dynamic outcomes associated with different dispersal behaviour. Integrodifference equations (IDEs), which are discrete in time and continuous in space, incorporate a dispersal kernel for the spatial distribution of mosquitoes, model populations in which growth and dispersal do not happen at the same time (Kot and Schaffer 1986;Hardin et al. 1990;Kot 1992;Neubert et al. 1995;Kot et al. 1996;Neubert and Parker 2004;Hsu and Zhao 2008;Zhou and Kot 2011;Reimer et al. 2016). Unlike the model approaches discussed in Table 1, IDEs allow a wide range of redistribution kernels for dispersal to be considered. This way of formulating time and space also provides a better approach for determining invasion speeds (compared to the reactiondiffusion PDE approach which can often underestimate patterns of spread; see Kendall 1965;Murray 2001).
The aim of this paper is to use optimal control theory to find the most efficient release strategy under different environmental conditions and dispersal behaviours. This is of critical importance as it can be used as a guideline for the number of sterile mosquitoes that need to be reared as well as aiding best practice for their release. Specifically, we develop a bioeconomic model with the corresponding cost function for the control  Joshi et al. 2007;Martinez et al. 2015) has considered optimal harvesting in an integrodifference population framework, but to the best of our knowledge no one so far has investigated the effectiveness of SIT through an IDE framework. In particular we will analyse the effect of spatial heterogeneity on the sterile release strategies. We find that an optimal control release strategy exists and efficiently suppresses the wild mosquito population. Our approach allows us to find the timing and the intensity of the control that needs to be applied. In Sect. 2 we introduce the integrodifference framework. We present the model in Sect. 3 and derive the characterization of the optimal control in Sect. 4. This characterization corresponds to the most cost-effective release strategy to suppress or eliminate the wild mosquito population. In Sect. 5 we use numerical simulations to illustrate the theory developed.

Integrodifference equations
In order to familiarise the reader with the integrodifference equation (IDE) framework (Kot and Schaffer 1986;Hardin et al. 1990;Kot 1992;Neubert et al. 1995;Kot et al. 1996;Neubert and Parker 2004;Hsu and Zhao 2008;Zhou and Kot 2011;Reimer et al. 2016), we give some basic notations and assumptions. Let N t be the population density at time t, f (N t ) the growth function where f (N t ) = N t g(N t ) and g(N t ) is the per capita growth. We model the population (in discrete time) with no movement as with initial condition N t 0 = N 0 for N 0 > 0. The spatial spread is described by a dispersal kernel, denoted by k(x, y), which is a probability density function (pdf) and gives the probability that an individual starting at point y, will settle at point x by the next time step. The number of individuals moving to location x is found by integrating the dispersal kernel k(x, y) over the domain of interest. Hence, we have k(x, y)dy ≤ 1 (as we are looking at a population on a finite domain), where is the spatial domain. We get the full dynamics of the population by combining the growth function with the dispersal kernel as follows:

The model
For the SIT control problem, we assume an unstructured mosquito population with no overlapping generations. Let W t denote the wild mosquito population at time t.
In the absence of sterile mosquitoes and (initially) ignoring any spatial variation, the governing difference equation for the wild mosquito population is where α is the number of matings per individual and β is the number of offspring produced per mating. γ (W t ) is the survival probability and we assume a Ricker-type nonlinearity, such that: The death function is approximated by the linear function d + K W t , where d is the density independent death rate, e −d is the density independent survival probability and K is related to the carrying capacity. Now define A = αβe −d (as the intrinsic population growth rate) and assume that no mating difficulties arise. Therefore, Eq. (3) takes the form: Suppose now that R t sterile mosquitoes are released at time t into the wild mosquito population. The total population size is then N t = W t + R t . The number of offspring produced by a wild mosquito over its life that would make it to adulthood without density-dependence, given that random mating with sterile mosquitoes occurs is AW t W t +R t . We assume that the life cycle happens in one time unit and that density dependence occurs only at the larval stage (Clements 1992;Lord 1998). This means that the survival probability depends only on the number of wild mosquitoes (and is independent of R t ). We assume that the number of released mosquitoes is R t = r t W * , where r t is the release ratio at time step t and W * is the equilibrium population; the release of sterile mosquitoes is proportional to the wildtype equilibrium population size. Under these assumptions the population dynamics in a non-spatial system are governed by the following equation: Next, we introduce spatial effects by incorporating the dispersal kernels and growth functions into our IDE framework. Assume that the wild mosquitoes have a dispersal kernel k(x, y) in a one dimensional domain . We assume a closed domain, for example an area that is surrounded by unfavourable conditions such that mosquitoes would not travel across. It is important to note that this might not always be the case and different spatial analyses will be needed for these different boundary conditions. To analyse the dynamics of the wild mosquito population under the influence of SIT release, we use: where f (W t (y), y) = AW t (y)e −K W t (y) and initial condition W 0 ( Table 2). The equilibrium population W * is calculated as

Optimal control formulation of integrodifference equation
In this section, we describe the optimal control framework for minimizing the cost of controlling the wildtype mosquito population. We propose to control a vector population over a time period representing the wild mosquito population and the sterile insect release ratio respectively at location x and time step t, where the initial distribution W 0 (x) is given.
The kernels are bounded and measurable such that | k(x, y)dy| ≤ 1, ∀x ∈ and 0 ≤ k(x, y) ≤ k 1 for (x, y) ∈ × , where k 1 is a constant. We assume that the We want to find the optimal strategy that suppresses the mosquito population and minimizes the cost of vector control via the release of sterile mosquitoes. Assume there is a linear cost associated with the wild mosquitoes (due to impact on human health, lost tourism, etc.) which we denote by m t . There is also a cost for producing and releasing sterile mosquitoes and we assume it to be a quadratic of the form,(n t r t W * + s t r 2 t (W * ) 2 ). The cost function is a nonlinear relationship in r*. This choice of cost function is based on the reality that sterile insect releases, through mating disruption, introduce an Allee effect into the wild population (Bonsall et al. 2010; also see Kirschner et al. 1997 for similar choice of quadratic costs for a HIV chemotherapy application). Using these costs, we define the objective functional: with set of bounded controls given by Our goal is to minimize the total cost of mosquito management by finding r * such that: In order to achieve this, we need to link the sensitivity ψ t (x) and adjoint system (append the difference equations (Eq. 6) to the objective functional to be minimized), to characterize the optimal control. By differentiating the state equation we show that the sensitivity satisfies: Proof of this result is given in the "Appendix". Using optimal control theory, we derive the adjoint system and the resulting optimal release strategy in the following theorem.
Theorem 1 Given an optimal control r * and corresponding state solution W * = W (r * ), there exist a weak solution λ ∈ (L ∞ ( )) T satisfying the adjoint system with the transversality condition and furthermore where X is the solution of and b and d are defined by: Proof of this Theorem is given in the "Appendix". We find the characterization of the optimal control by solving Eq. (13). Using the proof of Theorem 1 (see "Appendix"), if we analyse Eq. (32), we find its derivative is positive for all time. This implies that Eq. (32) has precisely one real root. If we obtain a negative real root, we set the control to zero. The solution is zero when C > A 1 indicating that the optimal strategy is to not release any sterile mosquitoes. This means that the burden imposed by the wild mosquito population (in terms of disease or nuisance biting) is so low that we do not need to instigate a sterile mosquito release programme. When the real root is positive, we release sterile insects at the release ratio dictated by the root.

Numerical results
We use a forward-backward sweep numerical method (e.g., Lenhart and Workman 2007) to find the optimal control r ∈ , where = (r 0 (x), . . . , r T (x)) ∈ (L ∞ ( )) T |0 ≤ r k (x) ≤ r max , t = 0, 1, . . . , T such that J (r * ) = min r ∈ J (r ). This iterative process starts by guessing the control (in our case we start with r = 0), and using the initial conditions of the state we solve the state equation forward in time. The next step is to use the new values of the state to solve the adjoint equation backwards in time. The values of the adjoint equation are used to solve the characterization of the control in order to get a new control value at each time step. To speed up the convergence, the value of the control is updated by taking a convex combination of the new and old control value. The new calculations of the optimal control and states are compared with results from the previous iterations. If the difference between iterations is less than 1%, the final control r is considered to be optimal.

Homogeneous environments
Here, we apply the optimal control theory to homogeneous environment. To replicate such environments we assume that the value of K (related to the carrying capacity) is the same across the whole domain . Throughout this section, we assume a uniformly distributed wild population at t = 0 with m = n = 1, s = 5 and consider = [0, 1], unless otherwise stated. We have chosen these values to correspond to a high cost of producing and releasing sterile insects compared to the cost associated with the wild mosquitoes. With this set of costs, we do not expect a complete elimination of the wild population. In the following, we investigate the Laplace kernel for the dispersal, which is described by: where the parameters associated with this kernel follow κ L = D L = 1. We have chosen an exponential form for the dispersal kernel, as these dispersal conditions have received the most support from the available data (Gratton and Zanden 2009;Estep et al. 2014). Model simulations show that in the absence of control, the wild mosquito population increases until it approaches the stable equilibrium at time step t = 5 (see Fig. 1). The distribution of the wild mosquito population changes through time. Initially it has a uniform distribution, and as a result of their redistribution and Fig. 1 The density of the wild mosquito population without control increases initially until it approaches the equilibrium at time step t = 5. A Laplace kernel is assumed for the dispersal in the domain = [0, 1] for A = 2 and K = 0.001 the chosen boundary conditions, mosquitoes tend to aggregate more in the centre of the domain as opposed to its boundaries.
When the control is introduced, the wild mosquito population decreases, but does not go to zero (see Fig. 2b). Similar to the behaviour without control, the population is aggregated in the centre of the domain and reduces towards the boundaries. For an effective SIT technique (capable of suppressing the wild mosquito population for the cost function and parameters given above), we need approximately 1.23:1 ratio of sterile to wild mosquitoes initially and releases should be applied more in the centre of the domain. The number of releases decreases through time (Fig. 2a) as the wild mosquito population size decreases to the point that they do not impose a high burden to the environment. This is the reason why the wild population size does not go to zero in Fig. 2b. Releasing sterile mosquitoes can reduce significantly the wild mosquito population, (by approximately 73% by time step t = 5) which could have a positive impact on disease management.

Heterogeneous environments
Very often, a habitat has different attractiveness to mosquitoes in different areas. This can be influenced by factors such as resource availability (e.g. food, mates, breeding sites) and predation. In this section, we analyse how landscape heterogeneity affects the dispersal and the optimal control of mosquitoes. We replicate this behaviour by varying the value of K (the strength of density feedbacks related to the carrying capacity). We consider situations (a) when the centre of the domain has more favourable Conditions (see Figs. 3, 4a), (b) when the boundaries of the domain have more favourable conditions (see Figs. 5,6a). We conclude that through time, the population in situation (a) aggregates more in the centre of the domain, as expected, because the conditions are more favourable there. This is true for populations with and without control. Without control, the wild mosquito population slightly decreases initially at the boundaries until it approaches the stable equilibrium. On the other hand, in the centre of the domain, the wild mosquito population increases until it approaches its equilibrium at t = 5. Once the control is introduced, we see from Fig. 4b that releasing sterile mosquitoes greatly suppresses the wild mosquito population. The optimal strategy in Fig. 2 a The number of sterile mosquitoes released over time in a homogeneous environment. b The density of the wild mosquito population under the optimal control strategy in the domain = [0, 1] for A = 2 and K = 0.001. A Laplace kernel is assumed for the dispersal. Releasing sterile mosquitoes does not eliminate the wild population, but it significantly reduces it Fig. 3 The density of the wild mosquito population increases at the centre of the domain. Wild mosquito population through time, with no control, when only the centre of the domain has favourable conditions. Laplace kernel is assumed Fig. 4 a The number of sterile mosquitoes released over time in the heterogeneous environment. b The density of the wild mosquito population under the optimal control strategy in the domain = [0, 1] for A = 2 and K = 0.001, when only the centre of the domain has favourable conditions. A Laplace kernel is assumed for the dispersal. Releasing sterile mosquitoes does not eliminate the wild population, but it significantly reduces it  The density of the wild mosquito population under the optimal control strategy in the domain = [0, 1] for A = 2 and K = 0.001, when only the boundaries of the domain have favourable conditions. A Laplace kernel is assumed for the dispersal. Releasing sterile mosquitoes does not eliminate the wild population, but it significantly reduces it this scenario is to release more sterile mosquitoes in the centre of the domain where they are concentrated. Figure 4a gives the optimal number of sterile mosquitoes that need to be released in order to control the wild mosquito population. We release substantially more mosquitoes at the centre of the domain than its boundaries at time step t = 1. This is because dispersing mosquitoes are moving more towards the centre and at this level the burden imposed by the wild mosquito population is so high that a large number of mosquitoes need to be released in order to suppress the population. Once the population is suppressed, after t = 2, the ratio of releases between boundaries and centre of domain is decreased. We get the opposite behaviour when the conditions are more favourable in the boundaries. Without control, from Fig. 5, we notice that the population increases initially and it approaches an equilibrium at time t = 5. As expected, the population size increases more in the boundaries, because mosquitoes will move more towards them as the conditions for mosquito reproduction are better. Using the control, we see from Fig. 6b that the population quickly decreases until it reaches a threshold where they do not impose a high burden. The optimal strategy in this case is to release more sterile insects in the boundaries initially (until t = 3) and after the wild population has reached a critical low threshold (e.g., where the effects are not very harmful), the ratio of release between boundaries and the centre of domain is decreased.

Discussion
Here, we formulated a novel mathematical model to understand the effects of releasing sterile mosquitoes into wild populations as well as the effect of spatial spread on mosquito population dynamics. The approach described here has not been used before in designing and optimizing sterile insect release strategies.
The model is described by an integrodifference equation, which are used to model populations where growth and dispersal do not happen at the same time. In our model we consider a homogeneous population consisting of wild and sterile mosquitoes with no overlapping generations. The growth function is based on the Ricker model and we assume releases proportional to wild population equilibrium. The spatial spread of the mosquitoes is described by the dispersal kernel k(x, y) in the integrodifference equation.
Using numerical simulations, we considered homogeneous environments with uniform carrying capacity, and heterogeneous environments with different carrying capacity in different areas, where the Laplace kernel describes mosquito movement. One significant finding is that, due to redistribution, applying the optimal control does not eradicate the wildtype mosquito but only substantially reduces population size. In practice, sterile mosquitoes may be released from aircraft resulting in releases that are approximately uniform. Our results highlight that this is not optimal and instead releasing more where the population densities are higher is more efficacious for vector control. Our model predicts a 73% suppression of the wild population which is close to observed field estimates, where > 80% suppression rates have been reported for A. aegypti control (Harris et al. 2012;Carvalho et al. 2015).
In heterogeneous environments, we consider situations where the centre of the domain has more favourable conditions or when the boundaries have more favourable conditions. In both cases the control significantly suppresses the mosquito population. Our results suggest (as expected) that more mosquitoes should be released where densities are higher. Another (expected) result from this model is that the optimal strategy is to release significantly more sterile mosquitoes at the beginning of the vector management control programme until the wild mosquito population is suppressed to a level that the imposed burden is so low.
Here we show that continuous releases predict that complete eradication of the wild population is not an optimal solution. Furthermore, integrating different control options (insecticide knockdowns, pulse or continuous SIT) to achieve cost-efficient control strategies needs more thorough investigation (Hackett and Bonsall 2018). Additionally, we argue that this discrete-time, continuous-space model approach is better than previous ones (described in Introduction) in finding the time and intensity of control, as it can include a variety of dispersal behaviour. This is crucial when modelling mosquito control. However, we emphasize the importance of determining a more accurate dispersal kernel that supports the field data, as the wave speed is very sensitive to the dispersal behaviour. Depending on the specific system parameterization, the results presented here are likely to be sensitive to the cost function form and the parameters used in this function. We have assumed a quadratic form, but other functions can be explored (see Khamis et al. 2018). This will change the objective functional J (r ) and the characterization of the optimal control. Importantly, our results emphasize that optimal control does not necessarily lead to population elimination. Varying the parameters (m, n, & s in Eq. 7) associated with the cost of the wild and sterile mosquitoes will most likely modify this outcome. If we lower the cost of producing and releasing sterile mosquitoes, we can achieve elimination of the wild population. In summary, we showed that the optimal control of the SIT model described by an integrodifference equation exists and that the control can significantly suppress the wild mosquito population.

By iteration we have
Subsequently we have W t (x)−W t (x) −→ ψ t (x) weakly in L 2 ( ). Using the pointwise convergence we can show that ψ satisfies

Proof of Theorem 1
Let r * be an optimal control and W (r * ) the corresponding state. For the arbitrary variation l with (r * + l) ∈ for very small positive , let W be the corresponding solution of the state Eq. (6). As the adjoint equation is linear in λ, there exists a weak solution λ which satisfies Eq. (11). To characterize the optimal solution, we take the directional derivative of the cost function J (r ) with respect to r in the direction l.
Since J (r * ) is a minimum value we have [m t ψ t (x)] dx by adding λ T ψ T (x) (which is zero due to the transversality condition) and using the adjoint equation: allows Eq. (27) to be written as: We can transform Eq. (32) as: Letting in Eq. (33) we obtain the following cubic equation: