Discrete-time stochastic modeling and optimization for reliability systems with retrial and cold standbys

1.   Introduction

West Nile virus (WNV) is maintained in a cycle between mosquitoes and hosts. Culex mosquitoes species are the predominant vector of WNV in North America ^[38,41], and birds are the most important host species for WNV amplification. Although many mammals, including humans, develop WNV disease, they do not typically develop sufficient viremia to transmit the infection ^[43].

West Nile virus was discovered in Uganda in 1937, but not first detected in North America until 1999, when an outbreak in New York City caused 62 human cases of WNV disease, 7 human deaths, and dramatic mortality in American crows ^[37]. The virus quickly spread throughout North America, and has since been reported in every contiguous state in the U.S. ^[37] and all but two Canadian provinces ^[23]. WNV is the leading cause of mosquito-borne disease in humans in the continental United States ^[10,13,14]. Since its introduction to the U.S., WNV is estimated to have caused 51,000 clinical human cases of WNV disease, 2,300 human deaths, and 7 million human infections ^[37]. In addition, the virus has been detected in over 300 species of birds and has caused substantial bird mortality ^[12]. Although many bird species have rebounded since the initial introduction of WNV to North America, a few species, including the American Crow, have suffered persistent reduced survival ^[20,37]. Species with limited range or that are otherwise endangered are especially vulnerable to the negative impacts of WNV ^[37]. An avian vaccine is available and was administered to wild populations of endangered California condors and rare island scrub jays ^[26]. However, since vaccination of wild populations is challenging and expensive, and since there are no vaccines or preventative medications available for WNV disease in humans ^[10,37], mosquito management remains an important tool for WNV control and prevention.

Local governments use mosquito management and integrated vector management plans to control mosquitoes and limit the spread of disease ^[4,22,36]. These mosquito control strategies combine cost benefit analysis and risk assessments in order to implement timely, effective, and environmentally responsible mosquito control ^[4,22,36]. Control measures considered in this paper include larvicide, which is applied to stagnant water to kill immature mosquitoes, and ultra-low-volume adulticide sprays, which target adult mosquitoes as they fly. Although research supports the benefit of managing mosquitoes to prevent disease, and previous work has provided many insights into principles of effective mosquito management, many questions remain open ^[4]. For example, the magnitude and duration of mosquito suppression required to prevent an epidemic and general principles for guiding effective and efficient treatment schedules are current topics of interest in mosquito management ^[4].

Mathematical and statistical modeling helps to elucidate mechanisms that drive WNV epidemics ^{[1,25,32,42,46]}. Several works investigate the impact of species diversity on the risk of an outbreak ^[1,25,32,42]. Variable host competence, including the potential dilution effect of dead-end hosts ^{[1,25,32,42,44]}, competition between bird hosts for resources ^[32], and biting preferences of mosquitoes ^[32,42,44], are a few of the ecological factors that can impact the risk and dynamics of a WNV epidemic. In relation to variable host competence, ^[25] presents a tool for estimating the basic reproduction number of WNV in any bird community using phylogenetic mixed effect models and eBird citizen science data. A statistical data analysis, which includes many of the above ecological mechanisms, supports heterogeneous host competence as a significant factor impacting the circulation of WNV within the bird community ^[42]. Interestingly, this study does not support the presence of a dilution effect. Environmental factors, such as land use, temperature-dependent parameters, and seasonal effects, can also impact the risk of a WNV epidemic ^{[32,42,44,46]}. A comparative study suggests the potential impact of temperature-dependence on the risk of a WNV epidemic outweighs that of ecological variability ^[44].

Additional mathematical analyses consider the control of WNV and mosquitoes. For example, much work supports the idea that a WNV epidemic may exhibit backward bifurcation ^{[1,2,6,15,50]}, which implies the robust management techniques are required to reduce the reproduction number of the disease below one. Several papers investigate the impact of mosquito control on the risk of infection by evaluating how control measures will likely impact the value of the reproduction number ^[7,15,47]. The previous works do not provide a dynamic model of the control action. Instead, results are obtained by relating control efforts to model parameters. These works support the idea that it is possible to terminate a WNV epidemic by reducing the size of the mosquito population. However, results also suggest that larvicide alone may not be sufficient for this purpose ^[15]. While analysis based on the model basic reproduction number widely supports the idea that the mosquito-to-bird ratio is an important determinant of the risk of a WNV epdidemic ^[44,47], uncertainty and variability in the model parameters make it difficult to determine a threshold value of concern. In ^[44] a mosquito-to-bird ratio greater than 100 indicates a high risk of WNV outbreak across a wide range of ecologically and environmentally feasible parameter values.

Finally, several works evaluate the impact of mosquito and disease control strategies ^{[2,6,21,31,50]}. In ^[21,50], a control acts to impulsively reduce the size of the mosquito population at periodic or state-dependent application times. Although it is possible to compute thresholds for the persistence of the disease in these models, due to their complexity, numerical simulations are needed to visualize the relationship between the control schedule and disease persistence. Alternately, some researchers have used optimal control frameworks to characterize effective and efficient control strategies ^[2,6,31]. In ^[2,6], continuous optimal control problems, where the effect of pesticide and personal protection can be continuously adjusted, are posed. Results suggest that larvicide may be more effective than adulticide when used in isolation ^[2], and that personal protection is less effective than mosquito control ^[6]. The former result seems to be in disagreement with ^[15], where it was suggested that adulticide is essential for disease control, however, this difference is likely due to the objective functional used in ^[2], which attempts to minimize, among other things, the size of the mosquito population. A discrete-time optimal control problem for a WNV epidemic is formulated in ^[31]. This problem is unique in that it considers, along with adulticide and pesticide, bird protection (e.g. vaccination), as a potential control action. Here the objective functional is designed to minimize mosquito larvae, infected birds, infected mosquitoes, and a quadratic control cost. The optimal control strategy is determined for a variety of parameterizations. In this study, optimal controls include high levels of bird protection, and larvicide is more heavily applied than adulticide. However, as in ^[2], the emphasis on larvicide is likely influenced by the fact that minimization of the larva population is an explicit goal of the objective functional.

In the current work, we use an impulsive differential equation model to study optimal treatment schedules for the purpose of mosquito control and WNV management. Our model is similar to that presented in ^[21,50], in that pesticide controls act to instantaneously alter the value of state variables. However, we allow adulticide and larvicide levels to decay continuously through time post application, so that pesticide decay rates can inform the optimal treatment timing. Indeed, in contrast to previous work ^[2,6,31], we optimize both treatment timing and application levels. In addition, we consider the problems of controlling the mosquito population and controlling the WNV disease burden separately, so we can better understand how different control schedules promote distinct objectives. In section 2 we formulate a WNV transmission model that describes the interaction between bird and mosquito populations (eggs, larvae, adults) and the dynamics for larvicide and adulticide. We derive the basic reproduction number of the infection using the next generation matrix method ^[45]. In section 3, we formulate two optimal control problems which seek to balance the cost of insecticide applications (both the timing and the application level) with (1) vector control: the benefit of reducing the number of mosquitoes, and (2) disease control: the benefit of reducing the disease burden. We reformulate the optimal control problems as nonlinear optimization problems, derive adjoint equations and establish optimality conditions. In section 4, numerical methods are discussed. In sections 5, results of simulations are presented. In section 6, discussions and conclusions are provided. Supplementary section 7 provides detailed descriptions of the model parameters.

2.   Model formulation & stability analysis

2.1.   Model formulation

The model describes the dynamics of WNV in the bird host population and female mosquito population. The bird population is divided into three classes according to their disease status. These classes are susceptible birds, $H_S (t)$, birds that are infectious with WNV, $H_{I}(t)$, and birds that have recovered and are no longer susceptible or infectious, $H_{R} (t)$. Hence, the total bird population is given as

The mosquito population is divided into three groups: eggs, aquatic larvae and pupae, and adults. These groups are further divided according to infection status. Since WNV can be transmitted vertically in the mosquito population, we have separate compartments for eggs laid by susceptible mosquitoes, $E_S(t)$, and eggs laid by infectious mosquitoes, $E_I(t)$. The aquatic population is further divided into susceptible larvae and pupae $L_S(t)$, and infected larvae and pupae $L_I(t)$. The adult stage group is divided into susceptible mosquitoes $V_{S}(t)$, exposed mosquitoes $V_{E} (t)$, and infectious mosquitoes $V_{I}(t)$. Hence, the total mosquito population is given as

where

and

Finally, the larvicide and adulticide levels are given as $U_L(t)$ and $U_A(t)$, respectively. Note that these levels do not match the absolute levels of pesticides in the environment or the biological impact of the pesticides as measured by percent inhibition of adult emergence or percent mortality in adults. Instead, they reflect the constant by which the maximal pesticide-induced larva and mosquito death rates should be scaled to match data on percent inhibition of adult emergence and percent mortality in adults, as described in supplementary section 7.

The model dynamics are described by the following system of differential equations.

where $d = \frac{m_L r_s q_S}{\mu_V}-\mu_L-m_L$, so that $C$ is the steady-state density of larvae.

The above system of equations reflects the following dynamic processes: Healthy susceptible birds acquire infection following contact with infected mosquitoes at a per capita rate $p_{MH} bV_I\frac {H_S }{ N_H}$, where $b$ is the per capita biting rate of mosquitoes on birds, and $p_{MH}$ is the transmission probability of West Nile virus per mosquito bite. Infectious birds die due to WNV disease at the per capita rate $d_H$ or recovered at the per capita rate $g_I$. We neglect host birth and death in this simple model. Hence this model is best suited to study a WNV epidemic occurring mid to late summer when most birds have finished laying eggs for the season.

Susceptible and exposed mosquitoes lay eggs at the per capita rate $r_S$ and infected mosquitoes lay eggs at the per capita rate $r_I$. Mosquito eggs mature at a per capita rate $m_E$. A proportion $q_S$ of the susceptible-laid eggs hatch into live larvae and a proportion $q_I$ of infectious-laid eggs hatch into live larvae. A fraction $\phi$ of larvae descending from infectious mothers are infected.

Larvae mature into adult mosquitoes at the per capita rate $m_L$. Larvae are subject to density-independent natural mortality at per capita rate $\mu_L$ and density-dependent mortality, with carrying capacity $C$. Larvae are subject to larvicide-induced mortality with a maximal per capita rate of $k_{m1}$.

Susceptible mosquitoes acquire the infection following contact with the infected birds at the per capita rate $\frac{b p_{HM} V_S H_I}{N_H}$. Adult mosquitoes have a natural per capita mortality rate $\mu_{V}$ regardless of their infection status. Exposed mosquitoes progress to the infectious class at the per capita rate $k_L$. Adult mosquitoes are subject to adulticide-induced mortality with a maximal per capita rate of $k_{m2}$. Adulticide and lavicide levels decay exponentially with rates $g_L$ and $g_A$, respectively.

The model variables and parameters are summarized in Table 1. In the supplementary section 7, we describe the model parameterization. Parameter values and ranges are summarized in Table 3. A model diagram is given in Figure 1.

2.2.   Stability of the disease free equilibrium (DFE)

The disease free equilibrium (DFE) of the virus-vector model (2.1) is given by

where

and $H_0$ is the initial density of hosts.

From the Routh Hurwitz criterion, we see that within the disease-free model, i.e., within the four-compartment model which is derived from the full model by neglecting the diseased variables ($H_I, E_I, L_I, V_E$), the DFE is locally asymptotically stable provided that

i.e., provided that $d > 0.$ That is, the eigenvalues of the Jacobian matrix of the four-compartment model evaluated at the DFE have negative real part in this case.

The stability of $\mathcal{E}_0$ within the full model can be established by calculating the reproduction number $\mathcal{R}_0$ using the next generation matrix method ^[45] on system (2.1). Taking $H_I, E_I, L_I, V_E, V_I$ as the infected compartments and adopting the notation in ^[45], the Jacobian matrices $F$ and $V$ which account for new infections and transfer between infected compartments, respectively, are defined as:

where $n_1 = d_H+g_I, \ n_2 = \mu_L+m_L+\frac{dL^*_S}{C} = \mu_L+m_L+d, \ n_3 = k_L+\mu_V.$ We define $\mathcal{R}_0$ as the spectral radius of the next generation matrix $FV^{-1}$. It is given by

where we have used (2.4) and (2.7) to eliminate $H^*_S$ and $V^*_S$ from the expression. $\mathcal{R}_0$ can be interpreted as the expected number of secondary infections, in both mosquitoes and humans, produced by an infectious individual during the infectious period in an entirely susceptible population. Note that since $V$ is a nonsingular matrix with the Z-sign pattern and $F$ is nonnegative, the disease-free steady state is locally asymptotically stable when $R_0 < 1$ and unstable when $R_0 > 1$ ^[45]. We summarize this result in the following theorem.

Theorem 2.1. The disease-free equilibrium $(H^*_S, 0, 0, E^*_S, 0, L^*_S, 0, V^*_S, 0, 0, 0, 0)$ of system (2.1) is neutrally stable when $\mathcal{R}_0 < 1$, and unstable when $\mathcal{R}_0 > 1$.

As parameterized (see supplementary section 7) our model has a basic reproductive number of $1.23$.

2.3.   The endemic equilibrium

We note that the model does not have an endemic steady state since new birds are not recruited into the population and vertical transmission in mosquitoes is insufficient to sustain the epidemic due to reduced fertility of infected mosquitoes. Indeed, suppose there exists a nonnegative endemic equilibrium $(H^i_S, H_I^i, H_R^i, E^i_S, E_I^i, L^i_S, L_I^i, V^i_S, V_E^i, V_I^i, 0, 0)$ setting all derivatives to zero, we see that $H_s'(t) = 0$ and $V^i_I\neq 0,$ yields $H^i_S\equiv0$, and thus $H^i_I\equiv0$ and $V^i_E\equiv0.$ Using the differential equations for $E_I$ and $V_I$ to solve for $E^i_I$ and $V^i_I$ in terms of $L^i_I$ and substituting the expression for $E_I^i$ into the differential equation for $L_I,$ we find

Note that (2.8) and $L_I^i, L_S^i\geq 0$ implies

Since $\mu_Vd = q_Sr_Sm_L-\mu_L\mu_V-m_L\mu_V > 0,$ nonnegativity of $L_S^i$ requires $\phi q_Ir_Im_L-\mu_L\mu_V-m_L\mu_V\geq 0$. Using the differential equations for $E_S$ and $V_S$ to solve for $E^i_S$ and $V^i_S$ in terms of $L^i_S$ and substituting the expression for $E_S^i$ and $L^i_I$ into the differential equation for $L_S,$ we find

From the previous discussion, we see that $q_Ir_I < q_Sr_S$ together with (2.10) implies $L_S^i < 0,$ a contradiction. On the other hand, $H_s'(t) = 0$ and $V^i_I = 0,$ immediately yields $E^i_I\equiv H^i_I\equiv V_E^i\equiv L_I^i\equiv 0$. Thus, as parameterized, there is no endemic steady state in any case. We summarize this result in the following theorem.

Theorem 2.2. If the only impact of WNV infection on the mosquito population is to alter the fertility of infected mosquitoes, the late-season model of WNV infection has a nonnegative endemic equilibrium if and only if WNV infection increases the fertility of infected mosquitoes.

3.   Optimal control problems

Mosquito and WNV controls are accomplished via a control schedule that is determined by $N$ treatment times, $T(i), \ i = 1\ldots, N$, spanning a fixed duration of time $T_f$, and the corresponding application levels of larvicide, $u_L(i), \ i = 1\ldots, N$, and adulticide, $u_A(i), \ i = 1\ldots, N$. Hence, the level of larvicide and adulticide in the environment is subject to discrete impulses:

The treatment times are discrete state variables. In order to select optimal treatment times, we introduce control variables $\tau(i)$, which give the time between treatment $i-1$ and treatment $i$, that is $\tau(i)$ is the $i^{th}$ waiting time. Hence,

and we have the terminal state constraint:

Thus, an impulse control for the system takes the form of a sequence of control triples $\left\{(u_{L}(i), u_{A}(i), \tau(i))\right\}_{i = 1}^N$ such that

where $t(i)$ is the minimal permissible value of the $i^{th}$ waiting time. In the sequel, $t(1) = 0$ and $t(i) = 1$ for $i = 2, \ldots, N$

We consider two objective functionals, P1 and P2, which are designed to promote mosquito control and disease control, respectively.

(P1: Vector Control) This objective functional seeks to balance the cost of insecticide applications with the benefit of reducing the number of mosquitoes. That is, we minimize the mosquitoes during the entire treatment period and the eggs at the final treatment time, while minimizing a quadratic control schedule cost.

We solve this impulse optimal control problem by first converting it to a nonlinear discrete optimization problem. For this, we introduce discrete state variables corresponding to the values of the continuous state variables just after the control application.

In addition, we introduce an artificial discrete state variable to record the value of the integral part of the objective functional,

and a discrete state variable to track the current time

After discretization, the objective functional becomes

To complete the discretization of the problem, we must define the transfer function, $G$, that determines the value of the state at the next time step in terms of the control and value of the state at the previous time step ^[8]. For convenience, we first define a function $g:\mathbb{R}^{13}\rightarrow \mathbb{R}^{13}$ representing the derivatives of the continuous variables. That is, the first twelve components of $g$ are given by the right-hand side of the model system (2.1) and the thirteenth component of $g$ is given by $c_V$ times the integrand in (3.8). So, for example, $g_1(x) = -p_{MH}bx_{10}\frac {x_1 }{ x_1+x_2+x_3}$ and $g_{13}(x) = c_V(x_8+x_9+x_{10})$. Furthermore, let us denote the solution of the initial value problem,

by $X(t; x_0).$ Then $G(x, u, \tau):\mathbb{R}^{14}\times \mathbb{R}^2\times\mathbb{R}\rightarrow \mathbb{R}^{14}$ is defined as follows: For $j = 1\ldots 13,$

where

and $h:\mathbb{R}^3\to\mathbb{R}^{14}$ is defined by

Notice that $G$ is a function of the state at the previous time and the current pesticide application levels and waiting time.

Now we begin to form the necessary conditions for an optimal control of the discrete problem: Choose $(u^*, \tau^*)\in\mathbb{R}^{3\times N}$ to minimize

subject to the following constraints:

for $i = 1\ldots N$

Notice we have $14(N+1)$ equality constraints on the state components plus $1$ equality constraint on the $14^{\rm{th}}$ state component at the final time, and $5N$ inequality constraints for the control variables. Also, recall that $t(1) = 0$ and $t(i) = 1$ for $i = 2, \ldots, N$. Letting

we define the $i^{th}$ vector of equality constraints according to

so that the equality constraints can be summarized as

Similarly, we define the vector of inequality constraints $q(z)$ according to

So that the inequality constraints can be summarized as $q(z)\leq 0$. Also, let $J(z)$ be the objective functional. Then, in the language of nonlinear programming, we seek to solve the following problem:

(P1' Vector Control): Choose $z^*$ to minimize $J(z)$ subject to $r(z) = 0$ and $q(z)\leq 0.$

Adopting notation from ^[8], we form the necessary conditions for this problem as follows.

Theorem 3.1. If $z^*$ is a solution of (P1), there exists a nonzero vector $\lambda = (-1, y(0), \ldots, y(N), y(N+1))\in{\mathbb{R}^{14(N+1)+2}}$ (where $y(i) \in \mathbb{R}^{14}$, for $i = 0\ldots N$ and $y(N+1)\in\mathbb{R}$) and a vector $\mu = (\mu_1, \ldots, \mu_N)\leq0$ (where $\mu_i \in \mathbb{R}^5$) such that if the Hamiltonian for the system is defined as

then

where the derivative of the Hamiltonian is

From the necessary conditions of Theorem 3.1, we derive equations for the adjoint variables, $y(i)\in\mathbb{R}^{14},$ $i = 0, \ldots, N,$ by considering the derivative of the Hamiltonian with respect to the state variables. In particular, differentiating with respect to $x(i)$ for $i = 0, \ldots, N-1,$

Recall that

so $\frac{\partial{G_{1:13}}}{\partial x{_{1:13}}(i)} = \frac{\partial X}{\partial x(i)},$ where

and

Differentiating with respect to $x(N)$, we find the adjoint variable at the final treatment:

where $K: = y(N+1).$ In fact,

Differentiating with respect to $u$ and $\tau$, we find the optimality conditions:

Now we consider an objective functional aimed at controlling disease. (P2: Disease Control) This objective functional seeks to balance the cost of insecticide applications with the benefit of limiting the presence of disease. That is, we seek to minimize the infected mosquitoes and hosts during the entire treatment period and infected eggs at the final treatment time, while minimizing a quadratic control schedule cost.

As for the vector control problem, we solve this impulse optimal control problem by first converting it to a discrete optimization problem. Since the process is very similar as for the first problem, we omit the details. However, the reader should note that the dynamics of $x_{13}$ differ between the two problems. Also, for the disease control problem, the constraint on the adjoint variable at the final treatment time is

4.   Numerical methods

We use the necessary conditions to identify candidate optimal controls. The necessary conditions are solved numerically using a simple forward-backward sweep algorithm ^[30]: $(1)$ The current control is used to update the state and adjoint variables. $(2)$ The state and adjoint variables are used to compute a new control from the necessary conditions, as described in Theorem 3.1. $(3)$ A convex combination of the new and current controls is used to update the control value, and the next iteration begins. The algorithm converges when the relative errors in the state, adjoint, and control variables are small. Specifically, the convergence criteria for each variable $z$ is

where $z$ represents the new value of the variable, $z^*$ represents the current value of the variable, $\delta\leq{10^{-7}}$ is a small constant controlling the size of the relative error, and $N$ is the number of components in the vector $z$. In addition, when computing the new control, we enforce an upper limit of $T_f$ on each waiting time. Although not formally included in the problem statement, this limit is implied by the terminal constraint $T_f = T(N).$ We use several standard initial guesses for the control schedule to initiate the forward backward sweep algorithm. For the vector control problem we begin with $\tau(i)\equiv 1,$ $u_A(i)\equiv u_L(i)\equiv 0.$ For the disease control problem we begin with $\tau(i) = \frac{T_f}{N-1}\ i = 2, \ldots, N,$ $\tau(1) = 0,$ and $u_A(i)\equiv u_L(i)\equiv 1.$

Since each objective functional is a continuous real-valued function of the control variables which are taken from a compact domain, an optimal solution is guaranteed to exist. Solutions of the control problems and necessary conditions, however, need not be unique. Hence when solving for the necessary conditions we are only identifying candidate optimal controls. For some objective functionals and parameterizations we find evidence of the existence of multiple solutions of the necessary conditions and local minima of the control problem. For example, the disease control problem with $c_E = c_I = 5000$, $c_A = 10$, $c_l = 1$, and $c_T = 0.05$ appears to have multiple local solutions as the value of the unknown constant adjoint variable, $y_{14}$, which determines the duration of the control period, $T_f$, increases past $0.92.$ For $y_{14} > 0.92$ and when initialized with our standard initial guess, the forward-backward sweep algorithm converges to a solution with a fairly regular control schedule and $\tau(1) = 2.41$. For $y_{14} < 0.92$ and when initialized with our standard initial guess, the forward-backward sweep algorithm converges to a solution with an irregular control schedule with $\tau(1)\approx 0$. However, a similar irregularly-timed solution can be found for $y_{14} = 0.925$ when the algorithm is initialized with the solution for $y_{14} = 0.92.$ Both solutions are shown in Figure 12. The irregularly-timed solution achieves a lower value of the objective functional than both the regularly timed solution obtained with the standard initial guess and the corresponding fixed-time solution. Hence the irregularly-timed solution is a good candidate for the optimal solution. Unfortunately, this preferred irregularly-timed solution seems difficult to find without a good initial guess. Since we are interested in exploring the potential benefits of optimally selecting the treatment times, as opposed to numerical optimization algorithms, in the sequel we only consider parameterizations of the problem for which the forward-backward sweep algorithm converges to a candidate optimal solution, in the sense that its limit outperforms the limit of the corresponding fixed-time algorithm. We note that the algorithm seems to more readily converge to such solutions for smaller $T_f$, larger values of $c_V, c_I$ and $c_E$, and smaller values of $c_T.$

Finally, before viewing the numerical simulations with control, it is worth noting that, since pesticides do not instantaneously eliminate mosquitoes and larvae, optimal controls will not apply any pesticide at the end of the treatment time.

Detailed descriptions of the model parameterization are in supplementary section 7. All simulations are performed in MATLAB.

5.   Simulations and results

5.1.   No control

Figures 2–3 show the egg, vector, and host densities without any control and after the introduction of a small number of infected mosquitoes into the population. We see that in the absence of control, the epidemic is predicted to burn out after approximately four months. At the end of the epidemic, the bird population is reduced to approximately $50\%$ of its original size. Interestingly, even at the epidemics height, only a small proportion of the birds and mosquitoes are infected. This may be due to the short lifespan of the mosquito and the speed with which the infection progresses within the bird hosts.

5.2.   Control

In the following subsections, we show results on the control of the mosquito population and WNV epidemic as determined by two objective functionals which are designed to achieve either mosquito control or disease control. In so doing, we seek to address the following questions: How does the control and its impact on the system vary with the objective type? Is it beneficial to optimize the treatment timing? How does the control and its impact on the population vary with the type of pesticide used?

To address the third question, we model two specific types of larvicide (S-methoprene briquets ^[40] and VectoBac ^[35]) and a general adulticide. Both larvicides are used to control mosquito larvae in stagnant water, such as road-side gullies, however, S-methroprene briquets are a slow-release larvicide treatment that is designed to last for several months, whereas the duration of effectiveness for VectoBac is much shorter. Hence, we are interested to see how the control schedule, cost, and effectiveness varies with the type of larvicide used. Our control also includes an adulticide. Although there are many types of mosquito adulticides, unlike larvicides they are similarly effective and short lived. Indeed, adulticides are typically distributed via ultra-low-volume (ULV) aerosol sprays which target flying mosquitoes so that the time the spray remains airborne (typically 15 minutes to one hour) determines the duration of the treatment. Although ULV sprays do not last very long, the adulticides that they contain are highly effective at killing mosquitoes. Since adulticides and larvicides target different mosquito life stages and have very different dynamics, we are interested to see differences in their deployment under each objective functional.

To address the second question we simulate in parallel to each optimal control problem the corresponding optimal control problem with fixed treatment times. From here on, we will refer to the problems defined in Section 3 as optimal-time control problems to contrast with the corresponding fixed-time control problems. We set the fixed-time schedule to apply the first treatment at the initial time, and apply subsequent treatments at equally-spaced intervals over the duration of the control period. In each case, the fixed-time optimal control objective functional is that of the original control problem less the quadratic control timing cost ($\sum\nolimits_{i = 1}^N \tau^2(i)$). When interpreting the control timing cost it is helpful to note that the timing cost is minimized subject to the constraint $\sum\nolimits_{i = 1}^N \tau(i) = T_f$, when $\tau(i)\equiv{\frac{T_f}{N}}$. That is, regularly timed control schedules are less expensive. Table 2 shows the values of the optimal-time objective functional and fixed-time objective functional evaluated at controls that satisfy the necessary conditions (i.e., candidate optimal controls) for each problem considered.

5.3.   Vector control (P1)

Figures 4–5 show the densities of eggs, larvae, and vectors under vector control with S-methroprene and VectoBac, respectively. We see that eggs, larvae, and adult mosquitoes decrease dramatically as a result of the control.

Figure 6 shows the larvicide and adulticide schedules with S-methroprene and VectoBac, respectively. These figures reveal several interesting control features. First note that both solutions involve controls with highly irregular timing wherein high levels of adulticide and larvicide are applied consecutively during the fist two days. Table 2 confirms significant gains are made by optimizing the control timing for this objective functional. Indeed, the control schedules that correspond to the optimal-time control algorithm achieve lower optimal-time objective functional values than their fixed-time counterparts. In fact, even the fixed-time objective value evaluated at the optimal-time control is lower than the fixed-time objective value evaluated at the fixed-time control. This means that even if we neglect to account for the quadratic timing cost, the optimal-time control is more efficient than the fixed-time control. After the second day, the control timing is more regular, however, the amount of pesticide applied continues to vary. For example, both solutions involve a third high dose application of larvicide, after which larvicide and adulticide application levels are significantly lower. In addition, Vectobac dosing has an unexpected oscillating pattern wherein high and low larvicide application levels alternate. Finally, note that control schedules with S-methroprene and VectoBac show very different long-term patterns. Whereas very little S-methroprene is applied after the third treatment, significant levels of VectoBac continue to be applied for the remainder of the control period. In addition, although the adulticide treatment schedule is similar between the two larvicides, the amount of adulticide applied is lower with S-methroprene. Looking at Table 2 we see that the vector control objective functional values are lower with S-methroprene than with VectoBac, indicating that the longer-lasting larvicide is more efficient according to this objective functional.

5.4.   Disease control (P2)

Figures 7–9 show the densities of eggs, larvae, vectors and hosts under disease control with S-methroprene and VectoBac, respectively. We see that these controls effectively prevent the spread of the infection and largely prevent mortality in the bird hosts. Interestingly, these control schedules also allow for the preservation and partial recovery of the mosquito population. Hence, our results support the idea that it is possible to prevent the spread of a WNV epidemic without eliminating the mosquito population. Figure 10 illustrates the larvicide and adulticide schedules with S-methroprene and VectoBac, respectively, under the disease control objective functional. We see that pesticide application levels under disease control are much lower than those under vector control. In addition, whereas larvicide levels are greater than adulticide levels under vector control, the opposite is true under disease control. In summary these results suggest that moderate pesticide application can control the spread of disease without eliminating the mosquito population and that adulticide is important for disease control.

Looking at Figure 10, we see that aside from a high initial treatment near $t = 0$ the control timing is nearly regular with almost equally spaced weighting times. This suggests that beyond the first waiting time, treatment times are selected primarily to minimize the quadratic timing cost, and hence there is limited benefit in optimizing the treatment timing beyond the first treatment. Looking at Table 2, we see that the objective functional values for disease control under optimal-time control are lower than those under fixed-time control, however, in each case, the difference between the two values is not large.

Finally, comparing the results for S-methroprene and VectoBac in Table 2, we see that VectoBac is only slightly less efficient than S-methroprene at controlling the spread of disease. This result differs from the the previous case, where S-methroprene was found to be much more efficient than VectoBac at controllong mosquitoes. Moreover, looking at Figure 8 we see that VectoBac has significantly less impact on the mosquito population. Hence our results suggest that short-lived pesticides may be preferred for disease management since they are highly effective for this purpose and have smaller environmental impacts.

In Figure 11, the model-predicted instantaneous percent reduction in adult emergence achieved by two larvicides over time under the simplifying assumptions discussed in section A.3 is given.

6.   Discussion & Conclusions

Here we present a West Nile virus transmission model that describes the interaction between bird and mosquito populations (eggs, larvae, adults) and the dynamics of larvicide and adulticide. Because our model does not include bird recruitment, it is best suited to describe mid- to late-season WNV dynamics. We parameterize the model to reflect warm conditions, conducive to mosquito reproduction and disease transmission, and derive the basic reproduction number of the infection ($R_0 = 1.23$) as the spectral radius of the next-generation matrix. We formulate two optimal control problems which seek to balance the cost of insecticide applications (both the timing and application level) with (1) the benefit of reducing the number of mosquitoes, and (2) the benefit of reducing the disease burden. Both problems are reformulated as nonlinear discrete optimization problems. Control schedules which satisfy the corresponding optimality conditions are identified numerically using MATLAB. Since necessary conditions need not be sufficient, we refer to these control schedules as "candidate optimal control schedules."

Numerical experiments address three questions: (1) How does the control and its impact on the system vary with the objective type? (2) Is it beneficial to optimize the treatment timing? (3) How does the control and its impact on the population vary with the type of pesticide used?

We find that the objective functional has a significant impact on the candidate optimal control schedules. For the vector control problem, the candidate optimal control is front-loaded with two high-level pesticide applications on successive days. The candidate optimal control schedule for the disease control problem exhibits more regular timing. Also, the level of larvicide in the environment under vector control is much higher than that under disease control. These differences allow the mosquito population to persist at reduced levels under disease control. This observation can have important implications for managers seeking to control WNV while minimizing environmental impacts.

The vector control problem also demonstrates that there can be substantial benefit in optimizing the treatment timing. In this case, we find that the value of the objective functional for the candidate optimally-timed control is considerably lower than that of the corresponding fixed-time control. In fact, the candidate optimally-timed control outperforms the corresponding fixed-time control even when we neglect to account for the timing cost. In contrast, for the disease control problem, we observe very minor benefit to optimizing the treatment timing. However, simulations suggest that it is important to optimize the initial treatment time, and that fast initial treatment is needed in this case (results not shown). Since the fixed-time control necessarily applies the initial treatment at $t = 0$, the benefit of optimizing the first treatment time is not evident when comparing the fixed and optimal treatment schedules for this problem.

Finally, when comparing solutions using long-lasting S-methroprene briquets to those using VectoBac, we see that applications levels for VectoBac are higher than those with S-methroprene. As a result, S-methroprene is more efficient than VectoBac at controlling the mosquito population. This increased efficiency is evidenced by the considerable difference in the vector control objective functional values for these two larvicides. Nonetheless, both larvicides achieve near complete elimination of the mosquito population under the vector control objective functional. In contrast, the disease control objective functional value with S-methroprene is only slightly lower than that with VectoBac, suggesting that these two larvicides are similarly efficient at controlling WNV. Since time-averaged larvicide levels in the environment appear to be lower with VectoBac, and since the mosquito population is larger under VectoBac, this result suggests that short-lived larvicides are sufficient for disease control, and may even be preferable due to their reduced environmental costs.

It is important to note that the management model developed here has many limitations. For example, since we do not consider bird recruitment, this model is not appropriate for studying WNV disease in the early spring and summer when considerable bird immigration and egg-laying take place. Similarly, the model does not incorporate spatial features or seasonality. Nevertheless, we believe our mosquito management model is helpful and reasonably realistic over short time periods and smaller spatial scales. Detailed characterizations of pesticide dynamics and optimal selection of pesticide application times, represent major contributions of this work, which have enabled us to derive the above-mentioned results. Future work will be aimed at extending the model to additional seasons and considering the impacts of spatial effects.

Acknowledgments

Agusto was supported by the Strategic Environmental Research and Development Program under grant RC-2639. Ding and Leander's work were supported by National Science Foundation (Award 1757493) and Ding's work was also supported by Middle Tennessee State University Non-Instructional Assignment Grant.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

A.   Supplementary: parameters explained

Here we describe how parameter estimates are obtained from available data. It is important to note that none of the data were collected for the purpose of parameterizing this model. Hence, a parameter's value/range represents a rough educated guess.

A.1.   Mosquito demographics

Culex mosquitoes are the primary vector of West Nile virus in North America ^[38,41]. The mosquito life-cycle consists of four stages, eggs, larvae, pupae, and adults. Mosquito eggs require water to hatch, and both the larva and pupa life stages are aquatic ^[11]. In our model, we combine larva and pupae into a single life stage, which we call the larva phase, so that the life cycle is compartmentalized into three phases: eggs, larvae, and adults.

When considering mosquito demographics, it is important to note that demographic parameters, including survival probabilities and maturation rates, vary widely with environmental conditions ^{[16,18,25,38,39,49]}. Since our model does not incorporate time-varying environmental conditions, we use parameter values that reflect favorable reproductive conditions for Culex mosquitoes.

In favorable conditions, Culex pipiens eggs hatch in about two days ^[11,18]. Hence, we set the hatch rate ($m_E$) to $\frac{1}{2}\ \mathrm{day}^{-1}$.

In warm conditions, Culex larvae mature in one or two weeks ^[11,16]. Hence, we set the maturation rate of larva ($m_L$) to $\frac{1}{7}\ \mathrm{day}^{-1}$.

Culex mosquitoes deposit eggs in rafts of about 100-200 eggs ^[18,41]. However, WNV may reduce the size of the first egg raft by $50\%$ ^[41]. Since the average length of a gonotrophic cycle for a Culex mosquito in a residential habitat is $5-10$ days (and depends on the species) ^[24], and since mosquitoes in the wild are unlikely to survive multiple gonotrophic cycles ^[24,41], we set the rate of egg-laying of uninfected mosquitoes ($r_S$) to $150/8\approx 18.75\ \frac{\mathrm{eggs}}{\mathrm{mosquito}\ \cdot{}\ \mathrm{day}}$ and the intrinsic rate of increase of infected mosquitoes ($r_I$) to $100/8\approx 12.5\ \frac{\mathrm{eggs}}{\mathrm{mosquito}\ \cdot{}\ \mathrm{day}}$.

WNV infection can additionally impact the viability of eggs. Based on data from ^[41], we set the probability of hatching for an egg laid by an infected mosquito to $43\%$ ($q_I = 0.43$), and we set the probability of hatching for an egg laid by a susceptible mosquito to $56\%$ ($q_S = 0.56$).

Larval death rates are influenced by density-dependent factors, such as cannibalism ^[38], and density-independent factors, such as extreme rainfall, and temperature ^[27,39,49]. Due to density-independent factors, it seems likely that larval death rates will be higher in the field than in the laboratory. Based on a study of Culex larva survival in a forested area ^[34], we set the density-independent larval death rate ($\mu_L$) to $0.16\ \mathrm{day}^{-1}$, so that the fraction of larvae surviving to adulthood in the absence of density-dependent effects is $\frac{m_L}{m_L+\mu_L}\approx 44\%$ in our model. Meanwhile, density-dependent effects are described via a larval carrying capacity, which determines the number of larva that a square meter of land area can support. This parameter is estimated as the product of the number of larval habitats per meter squared, the volume of an average larval habitat, the fraction of larval habitats that are occupied, and the number of larva per unit volume of occupied habitat. Estimates of the number of larval habitats per meter squared, the average volume of a larval habitat, and the fraction of larval habitats occupied were taken from ^[33], where the distribution of mosquitoes in low-land Hawaii was quantified. Specifically, following data from ^[33], $0.015$ habitats per meter squared represents an area rich in larval habitats, $0.8\ \mathrm{L}$ is the average volume of a larval habitat, and $0.007$ is the fraction of habitats occupied by Culex mosquitoes. The average number of larva per liter of occupied habitat is set to $50,$ in agreement with data from ^[9] on the average density of Culex pipiens mosquitoes in man-made breeding sites in Italy. This gives a larval carrying capacity ($C$) of $0.0042\ \frac{\mathrm{indv}}{\mathrm{m}^2}$. For comparison, the average density of larva within occupied larval habitats (such as puddles and buckets) in the Gambia was reported to be $22\ \frac{\mathrm{indv}}{\mathrm{m}^2}$ in ^[19] and the density of birds in low-land Hawaii is about $0.0065\ \frac{\mathrm{indv}}{\mathrm{m}^2}$ ^[33]. The larval carrying capacity will vary greatly with geographic region, weather, and land-use. For this reason we suggest a larval carrying capacity range of $0.001 - 0.01\ \frac{\mathrm{indv}}{\mathrm{m}^2}$.

Since the lifespan of an adult Culex pipiens mosquito is about $10.4$ days in the wild ^[24], we set the adult death rate ($\mu_V$) to $\frac{1}{10.4}\ \mathrm{day}^{-1}$. We assume that resources are plentiful so that adult mosquitoes are not subject to density-dependent death.

A.2.   Transmission, progression, and recovery

Parameters that determine WNV transmissibility vary widely between Culex species and with temperature ^[39]. For example, the probability of WNV transmission depends on environmental conditions, viral titre, and species ^[25,39], and may vary widely between individuals ^[41]. Based on the data from the previous sources, we set both the probability of transmission from host to mosquito ($p_{HM}$) and the the probability of transmission from mosquito to host ($p_{MH}$) to $0.50$. The disease progression rate ($k_L$), which is the reciprocal of the time between inoculation and the onset of infectiousness, is set to $\frac{1}{10}\ \mathrm{day}^{-1}$ based on data from ^[41]. The infection-induced death rate of hosts ($d_I$) is set to $\frac{1}{5} \mathrm{day}^{-1}$ based on data from ^[28]. The duration of viremia provides an estimate of the duration of infectiousness. Based on data from ^[28], we set the host recovery rate ($g_I$) to $\frac{1}{7}\ \mathrm{day}^{-1}$. We set the biting rate ($b$) to $\frac{1}{5}\ \mathrm{day}^{-1}$ based on data from ^[38,39]. Finally, some research suggests that vertical transmission of West Nile Virus can occur in mosquitoes. We set the probability that a larva deposited by an infectious mosquito is infected ($\phi$) to $2.8/1000$ ^[3].

A.3.   Larvicide and adulticide dynamics

We modeled two types of larvicide: long-lasting S-methoprene briquets ^[40] and VectoBac ^[35], both of which may be used to control mosquito larva in stagnant water, such as road-side gullies. The maximal rate of larvicide-induced death ($k_{m1}$) and the rate at which the larvicidal activity of S-methropene briquets decays ($g_L$) in our model are estimated from the percent reduction in adult emergence ^[40]. In order to relate these model parameters and data, we make several simplifying assumptions and approximations: (1) We suppose that the larvicide-induced, per-capita death rate ($k_{m1}U_L$) and hatch rate ($m_Eq_SE_S$) are constant over the interval in which adult emergence is measured. (2) We approximate the density of larvae over an interval as the steady-state density corresponding to the given parameters. (Note that the accuracy if this approximation improves as the value of $\mu_L+m_L$ increases.) (3) We suppose the steady-state concentration of larvae is well-approximated by that in the absence of density-dependent death. (Note this assumption holds when either the hatch rate or the density-dependent death rate is small.) Under the previous assumptions, the density of adults emerging over a time interval of length $T$ in the absence of larvicide is

and the density of adults emerging over the same interval of time in the presence of larvicide is

Hence, if $p_1$ and $p_2$ are the percent reductions in adult emergence achieved $D_1$ and $D_2$ days after larvicide application, respectively, and $p_M$ is the maximal reduction in adult emergence achieved at application then

Note that while $p_M$ is not necessary for model parameterization, we compute it for comparison with data. Table 3 shows the data used in these computations. Figure 11 shows the instantaneous percent reduction in adult emergence when the model is thus parameterized and under the above mentioned simplifying assumptions.

When considering the larvicide parameters, it is important to note that $U_L$ represents neither the absolute nor the normalized concentration of larvicide in the environment. Similarly, $U_L$ does not measure larvicidal effectivness, that is, percent reduction in adult emergence. Hence, $g_L$ is not comparable to the decay rate of a larvicide or the decay rate of the insecticidal activity of a larvicide.

Adulticides are distributed via ultra-low-volume (ULV) aerosol sprays which target flying mosquitoes ^[5]. Hence, the time the spray remains airborne determines the duration of the treatment. Adulticide droplets optimally sized to target flying mosquitoes are predicted to remain aloft for about 15 minutes to an hour in calm conditions. Hence, we set the settling rate ($g_A$) to $24\ \mathrm{day}^{-1}$ so that the average time to settle is $1\ \mathrm{hr}$. ULV sprays are highly effective at killing mosquitoes. To determine the maximum adulticide-induced death rate, we used mortality rate data for caged Ochlerotatus sollicitans exposed to a synthetic pyrethroid adulticide, Biomist 31:66, via ULV spray in the field. Specifically, 1 and 12 hours post-application, the mortality rate was about $89\%$ and $95\%$, respectively, in mosquitoes $91.4$ meters downwind of the spray site at $50\%$ the maximum label rate. The mortality rate was $4\%$ in the control group. In rough agreement with this data, we set the maximum adulticide-induced mortality rate ($k_{m2}$) to $110.52\ \mathrm{day}^{-1}.$ At this rate the adulticide is effective and fast acting: In the absence of other types of mortality, approximately $77\%$ and $90\%$ of mosquitoes initially exposed to the half-maximal dose will die within 1 and 12 hours, respectively, and $90\%$ is the asymptotic mortality due to exposure at the half-maximal dose. Specifically, $k_{m2}$ satisfies

It is important to note that adulticide parameters vary with atmospheric conditions, and the presence of vegetative and man-made shelters can reduce the effectiveness of adulticides ^[5].