Managing disease outbreaks: The importance of vector mobility and spatially heterogeneous control

Management strategies for control of vector-borne diseases, for example Zika or dengue, include using larvicide and/or adulticide, either through large-scale application by truck or plane or through door-to-door efforts that require obtaining permission to access private property and spray yards. The efficacy of the latter strategy is highly dependent on the compliance of local residents. Here we develop a model for vector-borne disease transmission between mosquitoes and humans in a neighborhood setting, considering a network of houses connected via nearest-neighbor mosquito movement. We incorporate large-scale application of adulticide via aerial spraying through a uniform increase in vector death rates in all sites, and door-to-door application of larval source reduction and adulticide through a decrease in vector emergence rates and an increase in vector death rates in compliant sites only, where control efficacies are directly connected to real-world experimentally measurable control parameters, application frequencies, and control costs. To develop mechanistic insight into the influence of vector motion and compliance clustering on disease controllability, we determine the basic reproduction number R0 for the system, provide analytic results for the extreme cases of no mosquito movement, infinite hopping rates, and utilize degenerate perturbation theory for the case of slow but non-zero hopping rates. We then determine the application frequencies required for each strategy (alone and combined) in order to reduce R0 to unity, along with the associated costs. Cost-optimal strategies are found to depend strongly on mosquito hopping rates, levels of door-to-door compliance, and spatial clustering of compliant houses, and can include aerial spray alone, door-to-door treatment alone, or a combination of both. The optimization scheme developed here provides a flexible tool for disease management planners which translates modeling results into actionable control advice adaptable to system-specific details.

where T N denotes the following N × N matrix: 1/22 Similarly, for reflecting boundary conditions, we have where T N and T N are denote the following N × N matrices: and

Equilibrium vector distribution
In order to define the next generation and second generation matrices, as well as solve the corresponding eigenvalue problems, we require the equilibrium solution to the mosquito population dynamics defined in Methods I of the main text, which can be written in vector-matrix form as: Here, N v and Λ denote N 2 -dimensional column vectors with components N v i and Λ v i , respectively, and µ N 2 denotes the N 2 × N 2 diagonal matrix with entries along the diagonal given by µ i . The equilibrium mosquito population distribution, whose value in patch i we denote by N ve i , is given by the stationary solution to Eq. (6). The vector of equilibrium mosquito population values, denoted by N ve , is found to be 3 Next generation matrix, second generation matrices, and eigenvalue problem derivation Throughout this section, let 0 k and 1 k denote the k × k zero and identity matrices respectively, for any positive integer k, and let 0 denote the N 2 dimensional zero vector.

2/22
Following [1], we define the 4N 2 dimensional new infectious vector F and the 4N 2 dimensional transfer vector V as and Here, we use [xy] as a short hand to denote a vector whose i th component is given by x i y i . In F and V, compartments 1 to N 2 refer respectively to exposed vectors in sites 1 through N 2 , compartments N 2 + 1 to 2N 2 refer respectively to exposed hosts in sites 1 through N 2 , compartments 2N 2 + 1 through 3N 2 refer respectively to infectious vectors in sites 1 through N 2 , and components 3N 2 + 1 through 4N 2 refer respectively to the infectious hosts in sites 1 through N 2 . We note that we have defined F to be the rate and which new infectious hosts and vectors are generated, rather than the more standard convention where F represents the rate at which newly exposed but not yet infectious hosts and vectors are generated. Our convention allows us to focus on the spatial spread of agents who are actually able to propagate the disease. In the subsection which follows our derivation of the second generation matrices and R 0 eigenvalue problem, we show that the choice between the two conventions is immaterial for determining the basic reproduction number. The matrices F and V are defined to be the Jacobian matrices of F and V taken about the disease-free equilibrium. We find the following: where the the 2N 2 dimensional square sub-matrices F 21 , V 11 , V 12 , V 22 are defined by and 3/22 Note that we use N ve N 2 to denote the N 2 × N 2 diagonal matrix with diagonal entries N ve i defined in Eq. (7). The inverse V −1 is given by The next generation matrix F V −1 is thus given by: where F 21 V −1 11 and −F 21 V −1 11 V 12 V −1 22 are given by and If E v , E h , I v , and I h exposed and infectious vectors and hosts are introduced as a perturbation to a completely susceptible disease-free equilibrium state, the distribution of newly generated exposed and infectious vectors will be given by the next generation matrix's action on the 4N 2 dimensional column matrix ( (where the superscript 'T' denotes the matrix transpose) as follows: The eigenvalues λ of the next generation matrix can be found by solving the following characteristic equation: Letting M denote the following matrix the above characteristic equation can be written more compactly: The corresponding eigenvectors of F V −1 are determined by the relation

4/22
Note that the eigenvalues of (F V −1 ) 2 are given by the set {λ 2 | λ an eigenvalue of F V −1 }, with corresponding eigenvectors determined by where The matrices M and M have identical spectra (this follows from the property that AB has the same spectrum of BA for any two square matrices A and B), and so we see that the eigenvalues λ of the next generation matrix can also be determined by solving We refer to M and M as the "second generation" matrices. The basic reproduction number, R 0 , is defined as the spectral radius of F V −1 . The next generation matrix is a non-negative matrix, so R 0 is, itself, an eigenvalue of F V −1 such that all other eigenvalues are no greater in magnitude [2]. Furthermore, to the eigenvalue R 0 , there corresponds eigenvector(s) with all non-negative components [2], and these will be the only biologically relevant eigenvectors of interest (there is no biological meaning to a negative or complex valued distribution of hosts or vectors). The non-negative eigenvectors corresponding to R 0 represent the worst case scenario initial distribution of infected hosts and vectors which produces the largest asymptotic infected growth rate under the disease dynamics linearized about the disease-free equilibrium; e.g. if one introduces a given total number of infected vectors and hosts into a disease-free equilibrium system according to some chosen distribution, choosing a distribution which corresponds to a non-negative eigenvector with eigenvalue R 0 will cause disease levels to grow as rapidly as possible (ignoring short-lived transient effects). Equations (24) and (25) show that these worst case scenario distributions will be comprised entirely of infectious (as opposed to exposed) vectors and hosts. We do not use the next generation matrix itself in our analysis, but instead utilize the second generation matrices to determine R 0 by finding the largest non-negative solutions to either Eq. (23) or Eq. (27), and we find the corresponding worst case scenario spatial distributions of infectious host and vectors by solving the following eigenvector equations:

Alternative more standard convention
Under the more standard convention, the new infection vector F * and the transfer vector V * are defined as and

5/22
The corresponding next generation matrix F * V * −1 can be shown to be where and We thus find the action of F * V * −1 on a distribution of infected hosts and vectors to be the following The eigenvalue problem for F * V * −1 is given by The eigenvalue problem for (F * V * −1 ) 2 is thus found to be where and Comparing the above expressions to Eqs. (22) and (26), we see that M * = M, and that M * and M have identical specta, and thus conclude that M * , M, M * , and M all all have identical spectra and have a largest eigenvalue R 2 0 . The second generation matrices under this more alternative more common convention indicate the spatial propagation of exposed but not infectious host and vectors, while our convention indicates the spatial propagation of infectious hosts and vectors.

.1 Isolated sites
For a neighborhood comprised of a single isolated site (e.g. N 2 = 1), the hopping rate ω and Laplacian matrix L are irrelevant, and our disease system reduces to a basic single-patch SEIR model. In this case, the characteristic equations in (23) and (27) are equivalent and trivial to solve, and we find where the single-site uncontrolled, non-compliant, and compliant basic reproduction numbers are defined, respectively, as

Homogeneous systems
For the special cases of no control, 100% compliance, and 100% non-compliance, all model parameters are equivalent at every site in the neighborhood. Specifically, the diagonal death rate matrix µ N 2 and emergence rate vector Λ are given by the following: No controls applied to system µ N 1 N 2 , Controlled system, 100% non-compliance µ C 1 N 2 , Controlled system, 100% compliance, No controls applied to system Λ N 1, Controlled system, 100% non-compliance Λ C 1, Controlled system, 100% compliance.
Multiplying both sides of Eq. (7) by the matrix (µ N 2 + ωL), the equilibrium vector population N ve for homogeneous systems is found to be The homogeneous system second generation matrices are given by 100% non-compliant system where R 00 , R 0N , and R 0C are the single-site basic reproduction numbers defined in Eq. (41). To find R 0 , consider first the case of no controls applied to the system, for which the characteristic equation (23) reduces to where R 00 is the single-site uncontrolled basic reproduction number defined in Eq. (41). The non-zero eigenvalues of F V −1 can be found by using the properties of determinants to manipulate Eq. (46) into and solving the following alternative characteristic equation for the eigenvalues z of the matrix K: where we define and where the eigenvalues z of K are related to the non-zero eigenvalues λ of F V −1 by The Laplacian is a symmetric positive semi-definite matrix [3], and so the matrix K is also symmetric and positive semi-definite [4]. Therefore, the eigenvalues z of K must be real and non-negative [4], and this fact together with Eq. (49) implies that the non-zero eigenvalues λ of F V −1 must be real and satisfy Thus, if F V −1 has ±R 00 for an eigenvalue, then its spectral radius will be given by R 00 . Equation , and the fact that all Laplacian matrices have zero for an eigenvalue [3] implies that Det [L] = 0. We thus conclude that Det [K] = 0 and that R 00 is the spectral radius of F V −1 for a homogeneous uncontrolled system. Analogous arguments can be made for homogeneous 100% non-compliant and compliant controlled systems, from which we obtain the following result.
Inputting Eqs. (45) and (51) into the eigenvector equations (28) and (29), and utilizing the fact that L1 = 0 for any Laplacian matrix, the normalized eigenvectors I h and I v of the second generation matrices corresponding to the eigenvalue R 2 0 , regardless of control or compliance, are found to be

Infinitely fast hopping
For the case of infinitely rapid mosquito hopping, we consider the limit ω → ∞. This is a singular limit due to the fact that the fact that the matrices (µ N 2 + ωL) −1 and (p v 1 N 2 + µ N 2 + ωL) −1 diverge and become ill-defined as ω → ∞, and so mathematical analysis must be handled with some care. Here, we utilize the methods of Tien et al. outlined in [5]: we consider the quantities µ 0 /ω and (µ 0 + p v )/ω to be small parameters, perform Laurent expansions of (µ N 2 + ωL) −1 and (p v 1 N 2 + µ N 2 + ωL) −1 about these small parameters, respectively, and retain the lowest order terms of the expansions (order negative one) as matrix expressions valid in the ω → ∞ limit. We give only the results of this process and use them to obtain an expression for R 0 in the ω → ∞ limit, and refer the reader to [5] for details on performing the Laurent expansions.
Following [5], we obtain the following in the infinitely fast hopping limit: where the expression g denotes the average of a site dependent quantity g over the entire neighborhood, and the product 11 T is the N 2 × N 2 matrix filled entirely with ones. Equations (7) and (53) yield the following equilibrium vector population: Similarly, the second generation matrices defined in Eqs. (22) and (26) reduce to the following: The above matrix is comprised of N 2 identical rows, so it has rank 1, and thus has only one non-zero eigenvalue (with algebraic multiplicity one) [4]. The matrix trace gives the sum of all eigenvalues for any square matrix, so the single non-zero eigenvalue of M (or equivalently M) is found by taking the trace of M, and R 0 will be given by the positive square root of the result. From Eq. (56), we find We thus have M = M = R 2 0 11 T /N 2 , and the normalized eigenvector solutions to Eqs. (28) and (29) are found to be

Infinitely slow hopping
In the infinitely slow hopping ω = 0 limit, all patches decouple from one another. This is not a singular limit, and we can find expressions for the equilibrium vector population and the second generation matrices by setting ω = 0 in Eqs. (7), (22), and (26): Here and throughout the rest of this paper, the subscript (0) refers to quantities evaluated for ω = 0. The second generation matrices are diagonal in the no-hopping limit, and the spectral radii are consequently trivial to find: The eigenspaces of M (0) and M (0) corresponding to R 2 0 are degenerate: any vector which is distributed entirely within any subset of the non-compliant sites will be an eigenvector with eigenvalue R 2 0 (degeneracy vanishes only for the special case of a spatially homogeneous system). Consequently, we know that the worst case distributions of infectious vectors and hosts will lie entirely within the non-compliant sites in the no-hopping limit, but we have no way of distinguishing which, if any, of all possible eigenvectors are more "important" or "correct" practically, in terms of biological meaning and control strategies. This ambiguity can be resolved to some extent by considering perturbations to the no-hopping case caused by small but non-zero hopping rates.

Perturbation theory for finitely slow hopping rates
For the case of slow but non-zero hopping, we analyze our system using degenerate perturbation theory. Here, we give the necessary definitions, results, and interpretations. The details of the perturbation formalism can be found at a utilitarian level in Ref. [6], or at a more rigorous mathematical level in Ref. [7]. Perturbation theory will provide accurate results for our system when the parameters ω/µ 0 and ω/(µ 0 + p v ) are much smaller than unity. Under this assumption, one can find simplified approximate expressions for N ve , M, and M by writing the inverse matrices in Eqs. (7), (22), and (26) as Taylor expansions in ω, and retaining only the zeroth and first order terms. We define the perturbations δN ve , δM, and δM through the following equations: where the zeroth order unperturbed terms are defined as the no-hopping limit expressions in Eqs. (59) and (60), and the perturbations are linear in ω. The perturbations are given explicitly by the following formulas:

10/22
where N ve N 2 (0) denotes the diagonal matrix with diagonal entries given by the components of N ve (0) , and δN ve N 2 denotes the diagonal matrix with diagonal entries given by the components of δN ve . Similarly, we define the perturbation to the basic reproduction number through the following decomposition where R 0(0) is the no-hopping basic reproduction number given in Eq. (61), and where δR 0 is linear in ω. The perturbed basic reproduction number squared is the spectral radius of the perturbed second generation matrices. It and the corresponding perturbed eigenvectors of M and M are defined by the following perturbed eigenvalue problems, where terms greater than first order in ω are ignored: where δI v and δI h are linear in ω. In the above equations, the vectors I h (0) and I v (0) are as-of-yet undetermined members of the degenerate eigenspaces of M (0) and M (0) , respectively, corresponding to the eigenvalue R 2 0(0) (e.g. vectors that are distributed in some unknown manner entirely within the non-compliant sites).
Assume that exactly J ≥ 1 sites are non-compliant, let {i 1 , i 2 , ..., i J } denote the indices of the J non-compliant sites, and let {e 1 , e 2 , ..., e N 2 } denote the standard orthonormal basis for R N 2 . RStandard degenerate perturbation analysis (see Refs. [6] and [7] for details) shows that the unperturbed vectors I h (0) and I v (0) , as well as the perturbations δR 0 , δI v , and δI h are ultimately determined by the spectrum of the J × J-dimensional matrix W defined by the following elements: In words, W is formed by the matrix elements of δM (multiplied by an overall scaling factor 1/2R 0N ) for which both indices correspond to non-compliant sites, and for our system, these elements turn out to be equivalent to the corresponding matrix elements of δM. Letting degC(n) and degN(n) denote the numbers of compliant and non-compliant nearest neighbors connected to a site n, respectively, we find the following explicit expressions for the components of W: If sites i j and i k are nearest neighbors 0, otherwise, where the constants κ and ξ are given by 11/22 and ξ = 1 The eigenvalues of W are proportional to the parameter κ, and κ is itself proportional to ω/µ 0 , with the constant of proportionality determined by the non-compliant single site basic reproduction number R 0N and the ratios of the non-compliant vector death rate µ N and extrinsic incubation period p v to the natural death rate µ 0 . The parameter ξ can vary between one and two depending on the strength of control in compliant and non-compliant sites. Note that the matrix −W is diagonally dominant, symmetric, and has all real non-negative diagonally entries, and so it is positive semi-definite and thus has all real non-negative eigenvalues [8]. Let z denote the smallest of these eigenvalues. Consequently, the matrix W will have all real non-positive eigenvalues, the largest of which will be −z. Thus, the largest eigenvalue of W is real and non-positive. Further, because W has all non-positive diagonal entries and all non-negative off-diagonal entries, the matrix W + s1 J will be non-negative given for any non-negative constant s ≥ max k {|W kk |}, where 1 J is the J × J identity matrix. Thus, by the Perron-Frobenius theorem for non-negative matrices [8], the matrix W + s1 J has a positive real eigenvalue λ P F (the value of which depends on s) for which all other eigenvalues are less than or equal in magnitude, and for which there corresponds an eigenvector x P F that can be taken to have all non-negative components. From these facts, we conclude that λ P F − s is the largest eigenvalue of W to which x P F is a corresponding non-negative eigenvector. Thus, W has a largest, real, non-positive eigenvalue with a corresponding non-negative eigenvector.
Applying standard techniques from degenerate perturbation theory, the perturbation δR 0 can be shown to be the largest eigenvalue of W. From this, we conclude immediately that δR 0 ≤ 0, and that δR 0 is always proportional to κ. The J-dimensional eigenvectors α defined by the relation can be taken to be non-negative and normalized. Standard degenerate perturbation analysis shows that the components of α determine I h (0) , I v (0) , δI h , and δI v as follows: The first double summations in δI h and δI v are the perturbations to I h (0) and I v I v (0) . The matrix element perturbations δM ij and δM ij are non-zero only if i = j or if site i and site j are connected, and so the eigenvector perturbations δI h and δI v effectively transport the unperturbed eigenvectors from connected blocks of non-compliant sites out into the corresponding bordering compliant sites. It should be noted that, due to differences in normalization schemes, the analogues of the first double summations in Eqs. (74) and (75) are not present in many of the perturbation analyses found throughout the literature. Physics applications, for example, most often normalize vectors by setting their absolute magnitudes to unity, and for our system, this more typical normalization scheme would force the perturbations δI h and δI v to be orthogonal to I h (0) and I v (0) , thus eliminating the first double summations. The absolute magnitude normalization, however, does not have a clear biological meaning for our system. The alternative normalization scheme used in this paper has clear biological meaning in allowing I h and I v to represent probability distributions. 6 Relating controlled model parameters to real-world control strategies

Control strength
In our model, door-to-door adulticide and larval source reduction uniformly increase death rates and reduce emergence rates, respectively, in compliant sites only, while area-wide adulticide spray uniformly increases death rates in all sites. To door-to-door adulticide and aerial spray adulticide, we associate the control strengths ρ D ∈ [0, 1] and ρ A ∈ [0, 1], respectively, and to door-to-door larval source reduction, we associate the control strength σ D ∈ [0, 1]. A management strategy's control strength is dependent upon both the efficacy of the strategy as well as the frequency of application, and is defined as the percent reduction in equilibrium vector population that would result from that strategy's use in an isolated unconnected neighborhood site. By this definition, for a single type of adulticide used alone at control strength ρ i , i ∈ {D, A}, the increased vector death rate is given by µ 0 /(1 − ρ i ), and for larval control used alone at control strength σ D , the decreased emergence rate is given by Λ 0 (1 − σ D ), where µ 0 and Λ 0 are the natural death and emergence rates, respectively. Defining control strength in this way allows one to associate changes in model parameters, and therefore predictions for outbreak potential via R 0 , with notions of control effort; for high control effort, one expects large population reductions and control strengths close to 1, and for low control effort, one expects small population reductions and control strengths close to 0. The results in Ref. [9] provide explicit expressions for the aerial spray adulticide control strength ρ A (f A ) and the door-to-door adulticide and larval control strengths ρ D (f D ) and σ D (f D ) written as functions of the aerial and door-to-door application frequencies f A and f D , and are considered to be implicit functions of real-world measurable control properties like aerial spray percent knockdown, residual barrier spray efficacy time, and larval source recovery time, for example. The precise mathematical form of these expressions are given in the subsections below. Generally, control strengths are increasing functions of their application frequencies and approach zero as application frequencies approach zero. Numerical values for real-world control parameters appearing in control strength expressions are chosen based on expert opinion (Kevin Caillouet, personal communication) and the fitting results in Ref. [9], and are given in Table A. Figure A shows plots of ρ A , ρ D , and σ D as functions of application frequency under our assumed model and control parameters. Following the results in Ref. [9], the individual control strengths ρ A (f A ), ρ D (f D ) and σ D (f D ) combine 13/22 to give the following death and emergence rates in compliant and non-compliant sites:

Adulticide aerial spray
Area-wide spraying consists of an adulticide such as Naled, Malithion, or a pyrethoid which is dispersed throughout the neighborhood at ultra-low volumes by truck or airplane [10]. The spray remains suspended in the air for minutes to at most an hour, and kills adult mosquitoes upon contact [10]. We model adulticide aerial spray as an instantaneous fractional knockdown in mosquito population. Specifically, when adulticide aerial spray is applied to an isolated neighborhood site, the mosquito population is instantaneously reduced by the percent knockdown ρ. For repeated impulses applied at frequency f A , Ref. [9] shows the associated control strength ρ A (f A ) in Eq. (76) to be 14/22

Door-to-door adulticide
Residual barrier spray consists of a long-lasting adulticide such as bifenthrin or lambda-cyhalothrin, which is applied by individuals (with backpack mounted or hand sprayers, for example) to vegetation and other potential mosquito landing sites [10]. Adult mosquitoes who land on treated surfaces are quickly killed, and killing efficacy is typically retained for days to weeks as the insecticide evaporates away [10]. We model residual barrier spray as an instantaneous increase in vector death rate, and a subsequent recovery to the the natural value. Specifically, when door-to-door adulticide is applied to an isolated neighborhood site, the vector death rate is instantaneously increased by a fraction γ, and the increase decays exponentially at rate η (meaning that 1/η is the barrier spray efficacy time). For repeated impulses applied at frequency f D , Ref. [9] shows the associated control strength ρ D (f D ) in Eq. (76) to be where we define and where s > 0 and Γ denotes the doubly incomplete gamma function:

Door-to-door larval source reduction
Larval source reduction consists of employees identifying and eliminating receptacles for standing water typically found in residential yards (such as tires or empty buckets) which can serve as habitats for mosquito larvae. We model larval source reduction as an instantaneous decrease in vector emergence rate, and subsequent recovery to the natural value. Specifically, when door-to-door larval source reduction is applied to an isolated neighborhood site, the vector emergence rate is instantaneously decreased by a fraction σ, and the emergence rate recovers exponentially at rate ν (meaning that 1/ν is the carrying capacity recovery time). For repeated impulses applied at frequency f D , Ref. [9] shows the associated control strength σ D (f D ) in Eq. (77) to be Here, we present optimized control results for finite hopping rates analogous to those given in the main text, but for dispersed and clustered 20% compliance distributions under both periodic and reflecting boundary conditions. The plots are given given in Fig. B. At 20% compliance, for the clustered distribution in Fig. Bf under either boundary condition, optimal control action calls for combined strategies for hopping rates ω ∈ (6.20µ 0 , 10.5µ 0 ). For hopping rates below and above this interval, optimal control action calls for aerial spray only and door-to-door control only, respectively. The corresponding intervals for the dispersed 20% compliance distribution are given by (1.30µ 0 , 1.90µ 0 ) under periodic boundary conditions and (1.70µ 0 , 2.30µ 0 ) under reflecting boundary conditions. For each boundary condition and compliance configuration pictured in Fig. Bf, there exists a hopping rate interval where system is controllable under door-to-door control alone, but is more expensive to control with door-to-door control than with aerial spray. For hopping rates below this interval (excluding ω = 0), the system is uncontrollable with door-to-door control alone, and door-to-door only control is optimally applied as frequently as possible (once per day). Above this interval, the system is controllable under door-to-door only and is cheaper to control than with aerial spray only. The intervals are given by (6. . Boundary conditions are taken to be periodic. For each compliance configuration and control strategy, we simulate the full disease dynamics with initial conditions given by the susceptible disease-free equilibrium state plus a single infectious vector introduced into a non-compliant site. For each simulation, we choose ω equal to either 10µ 0 , 5µ 0 , µ 0 , 0.5µ 0 , or 0.1µ 0 , select an application frequency between 0 and 1 day −1 , compute R 0 , and find the cumulative total number of infectious hosts generated over the entire course of the epidemic. Results are displayed as curves depicting the total epidemic size as a function of R 0 . Each curve represents a particular hopping rate, and arrows along the curves indicate directions of increasing application frequency such that the leftmost endpoint of each curve corresponds to the maximum application frequency 1 day −1 . When this point lies to the left of the dotted R 0 = 1 line, the system is considered controllable (in the sense defined in the main text) under the corresponding door-to-door only or aerial spray only strategy. Figures C and D show that for a given set of system parameters and initial conditions, both R 0 and total epidemic size decrease with increasing application frequency. In other words, when controlling a system using either aerial spray only or door-to-door only strategies, total epidemic size decreases monotonically with the degree to which R 0 is reduced. For systems which are controllable, Figs. C and D show that reducing R 0 to or below one ensures the total epidemic to be small and contained.
For all controllable combinations of hopping rate and compliance distribution considered in Fig. C, with one exception, the combined cost-optimal R 0 = 1 controls as calculated in the main text call for either door-to-door only strategies or aerial spray only strategies. For these cases, the total size relations at R 0 = 1 under the cost-optimal strategy are pictured in either to the right of the R 0 = 1 line), the cost-optimal control is an aerial spray only strategy, and the corresponding R = 1 total size relations are given in Fig. D. The exceptional case is the 20% clustered compliance distribution with hopping rate ω = 10µ 0 (the teal curve in Fig. Cd). Here, the system is controllable under door-to-door only, and the cost-optimal control calls for a combined door-to-door and aerial spray strategy. The cost optimal total size at R 0 = 1 is found to be 11.9 hosts, which is comparable to the door-to-door and aerial spray only values of 15.0 and 16.2 hosts, respectively. We thus conclude that cost-optimal strategies for reducing R 0 to one coincide with strategies which reduce outbreaks to small sizes. It is important to note that cost-optimal controls for reducing R 0 to one are not not equivalent to optimal controls for reducing the total epidemic size. Specifically, combined aerial spray and door-to-door strategies which bring R 0 to one while minimizing the total epidemic size are unlikely to coincidentally minimize the cost of control. Additionally, cost-optimal controls for reducing the total epidemic size (rather than for reducing R 0 ) will be time-dependent, and finding them will require more sophisticated optimal control theory tools.  Fig. 7f of the main text. Teal, blue, red, green, and purple curves correspond to the hopping rates ω = 10µ 0 , 5µ 0 , µ 0 , 0.5µ 0 , and 0.1µ 0 , respectively, and the dotted vertical line indicates R 0 = 1. Each curve shows the total epidemic size as a function of R 0 for a fixed initial condition and various door-to-door application frequencies between 0 and 1 (day) −1 . Black arrows along curves indicate directions of increasing application frequency. In a given plot, the rightmost endpoints of all curves converge to the uncontrolled R 0 value 2.87 (representing an application frequency of 0), and the leftmost endpoints marked with large dots represent the maximum application frequency 1 (day) −1 . All figures show that the total outbreak size will be small whenever R 0 ≤ 1 . Aerial spray only control: Relation between the total epidemic size of infectious hosts and the basic reproduction number. Teal, blue, red, green, and purple curves correspond to the hopping rates ω = 10µ 0 , 5µ 0 , µ 0 , 0.5µ 0 , and 0.1µ 0 , respectively. Each curve indicates the total epidemic size as a function of R 0 for a fixed initial condition and various aerial application frequencies between 0 and 1 (day) −1 . Black arrows along curves indicate directions of increasing application frequency, and the dotted vertical line indicates R 0 = 1. The rightmost endpoints of converge to the uncontrolled R 0 value 2.87 (representing an application frequency of 0), and the leftmost endpoints marked with large dots represent the maximum application frequency 1 (day) −1 . All curves nearly completely overlap one another, indicating that the total size relation is insensitive to hopping rate for aerial spray only control strategies.