MODELING THE SPREAD OF BED BUG INFESTATION AND OPTIMAL RESOURCE ALLOCATION FOR DISINFESTATION

. A patch-structured multigroup-like SIS epidemiological model is proposed to study the spread of the common bed bug infestation. It is shown that the model exhibits global threshold dynamics with the basic reproduction number as the threshold parameter. Costs associated with the disinfestation process are incorporated into setting up the optimization problems. Procedures are proposed and simulated for ﬁnding optimal resource allocation strategies to achieve the infestation free state. Our analysis and simulations provide useful insights on how to eﬃciently distribute the available exterminators among the infested patches for optimal disinfestation management.

1. Introduction. Mathematical models are practically and conceptually valuable to address questions and challenges raised by complex demographic processes in ecology [22]. In urban systems, individuals travel daily between distinct spatial locations such as urban centers, cities, or local communities. The spread of societal pests such as the common bed bug, Cimex lectularius L., is largely influenced by the mobility of individuals.
The common bed bug is ranked as one of the most challenging of all insect pests to control due to the cryptic behavior and the insecticide resistance they have developed in recent years [11]. Since the 1990s, bed bugs have been considered widespread societal pests [10]. The number of infestations of the common bed bug and the risk of exposure through normal daily life have undergone a major resurgence worldwide and particularly in developed nations [2,16,23]. This has raised tremendous concerns and the challenge of pest management subject to a financially optimal criterion is of great interest to public health departments [8,26].
The majority of the bed bug infestations are associated with unsanitary living conditions and severe crowding, there are numerous evidences that infestations can occur in single-family dwellings, apartments, hotels [15], public transports such as trains [6], schools [26], furnished temporary housing including corporate apartments and serviced apartments [6] and hospitals [8]. It has been documented that cities with a high concentration of migratory factory workers experience frequent bed bug infestation [33]. The transmission process of the bed bugs occurs primarily when people visit infested sites or move the furniture and furnishings from infested sites [7]. Additionally, bed bugs can hitchhike on items such as clothing, backpacks and luggage [26].
Successful bed bug infestation management is extremely difficult due to the following reasons: (1) increased worldwide travel, which increases the risk of having multiple infestation sources [32]; (2) the cryptic behavior of the insects, which may allow them to remain undetected for several months until the infestation occurs on a large scale [14,25]; (3) its pesticide resistance, which can allow a reduction but not eradication in bed bug populations [18]. The integrated pest management method that combines the nonchemical means of control (such as heat) with the use of approved insecticides has been proven to be effective in the treatment of bed bug infestations [18,21].
Managing a bed bug infestation can be expensive in the costs of labor and chemicals. For instance, from 2001 to 2004, Australia experienced a resurgence in bed bug infestations involving both the common (Cimex lectularius L.) and Tropical (Cimex hemipterus) species. It was found that all Australian mainland states had experienced an exponential rise in bed bug infestations since 2001, with an overall national increase of 4.5%. A conservative estimate of the economic impact of the resurgence was 100 million Australian dollars [12]. A survey in New York during 2010 revealed that the average bed bug disinfestation cost was $1, 310, while across the United States, the average costs were around $800 to $1, 200 for the disinfestation treatment of a single unit [8].
The epidemiological approach has been employed to model many biological and ecological phenomena such as the viral propagation of memes [34]. There also has been some work on pest control management within a single patch (see, for example, [13,35]). In both [13] and [35], impulsive biological control strategies were proposed for pest control management purposes. However, no patch-structure or control related costs were considered in [13,35].
The purpose of this work is to understand the infestation dynamics of bed bugs and to explore optimal management strategies via a mathematical modeling approach. More precisely, we propose a patch-structured multigroup-like SIS system to model the processes of transmission and extermination of the common bed bugs, and we set up optimization problems to find the optimal resource allocation strategies to efficiently distribute the available exterminators among the infested patches for optimal disinfestation management.
The rest of the paper is organized as follows. We formulate our model and establish the global threshold dynamics in Section 2, with the proof of our main theorem postponed to the Appendix. In Section 3, we set up the optimization problems and propose procedures to find the practical optimal resource allocation strategies. In Section 4 the case of two patches is considered, basic reproduction number is computed and the optimization problems have been set up as an example and for illustration purposes. We summarize and discuss our findings in Section 5.

2.
Model formulation and analysis. Assume that there are n patches under consideration in a spatially homogeneous environment. Here a patch can be as large as a continent or as small as a local community. In patch i (i = 1, . . . , n), the total population N i is divided into two compartments: the susceptibles S i and the infesteds I i (thus N i = S i + I i ). Here by population we mean the number of units, which can be a house, an apartment unit, a dormitory or a hotel. At time t, a susceptible unit is a unit that has no bed bugs but can be infested. An infested unit is a unit which has bed bugs. It is assumed that if a unit is partially infested (for example only one room of a house is infested), then the whole unit will be considered to be an infested unit. This is due to the bed bugs' active movement through wall voids such as those used for wires or pipes allowing for a greater potential for other rooms to be affected by the bed bugs [26,27]. Since the transmission of the bed bug infestations primarily occur when people's daily activities involve visiting infested sites, it is reasonable to assume that the transmission occurs not only within a patch but also among patches. We use β ij to denote the transmission rate of the bed bug infestation from patch j to patch i (j = i) and β ii the transmission rate of the bed bug infestation within patch i. It has been documented [12] that in the absence of a thorough treatment, the bed bug infestation increases exponentially. Therefore, β ij s can be assumed to be positive constants. The incidence rate, the average number of new infested cases per unit time, in patch i due to the effective contacts from patch j is in the form of standard incidence β ij S i I j /N j . Thus, the overall incidence rate is the summation of contributions from all patches. The constant per capita extermination rate of the infested units in patch i is denoted by γ i . Thus, 1/γ i is the duration of the infestation period with an exponential waiting time in patch i. Once an infested unit becomes infestation free, it goes back to the susceptible compartment as it can be infested again. Thus, a repeat infestation right after an extermination treatment is possible [9].
The above assumptions lead to the following susceptible-infested-susceptible (SIS) model described by a system of 2n ordinary differential equations with initial conditions S i (0) > 0, I i (0) ≥ 0, and n j=1 I j (0) > 0. A schematic flow chart of the model with n = 2 is shown in Figure 1. It is convenient to work with proportions in all compartments given that N i (t) = N i (0) for i = 1, . . . , n. Rescaling the variables in system (1) as follows yields (3) Here in system (3), N i is rescaled to 1 for i = 1, . . . , n.
The non-negative transmission matrix of the bed bugs, F = (β ij ), can be further assumed to be irreducible. This means that the n patches cannot be separated into two groups such that there is no transmission of the bed bugs between these two groups. Mathematically, matrix F is irreducible if for every nonempty, proper subset J of the set N = {1, . . . , n}, there is an i ∈ J and j ∈ N \ J such that β ij = 0 [28].
Since the population size in each patch is fixed, system (3) can indeed be totally decoupled into two subsystems: the susceptible subsystem and the infested subsystem. The infested subsystem obtained by replacing S i with 1 − I i is of the following form where I = (I 1 , I 2 , . . . , I n ) T . Thanks to the results developed in the theory of asymptotically autonomous systems [3], it suffices to analyze the infested subsystem only.
Theorem 2.1. Consider system (4). If R 0 ≤ 1, then the IFE I 0 is globally asymptotically stable and if R 0 > 1, then I 0 is unstable and there exists a unique positive equilibrium I * , which is globally asymptotically stable.
3. Optimal resource allocation for disinfestation management. In this section we first identify the underlying parameters involved in the definition of the basic reproduction number R 0 in (5) and then explore optimal resource allocation strategies for disinfestation management by associating the disinfestation process with a cost function. To demonstrate our analysis, numerical simulations are carried out for an example of 2 patches.
Note that contacts do not necessarily transmit bed bugs [8]. For each contact made between two units, of which one is infested and the other is susceptible, there is a mean probability that the bed bugs can be transmitted. This probability is often called "transmissibility" [1]. We denote by p ij the transmissibility from patch j to patch i. Another factor that can affect β ij is the average number of contacts per unit in patch i with units in patch j per unit time which we denote by c ij . This can be interpreted as the average number of visits or travels from patch j to patch i per unit time. Thus, the transmission rate is The other parameters involved in the basic reproduction number are the extermination rates, γ i for i = 1, 2, . . . , n. The extermination rate of bed bugs in each patch depends on many factors such as the number of pest control companies that offer disinfestation services, methods of extermination, civilian knowledge of bed bugs, and media coverage. We choose the first two key factors to define the extermination rate γ i in patch i. In practice, each exterminator has several approved extermination methods, each of which takes a certain amount of time to fulfill the extermination process. Here an exterminator refers to an extermination unit which is equipped with required equipment and one or more licensed individuals who can perform the disinfestation procedure. Thus the extermination rate in patch i is given by where M ki is the number of exterminators in patch i using method k, T ki is the time taken by method k in patch i to exterminate the bed bugs, and the index m i represents the total number of extermination methods applied in patch i. Let C ki be the average cost of method k applied in patch i, then the total cost of the disinfestation of bed bugs among n patches is given by In practice, the extermination methods applied by exterminators to treat a typical bed bug infestation can be characterized into two types: (1) preliminary method and (2) advanced method. The preliminary methods refer to those at low cost but require a longer treatment time (relative to advanced methods), while the advanced methods refer to those at high cost but require a shorter treatment time. Examples of the preliminary method include the use of trained canines, regular and thorough inspections by residents of a site or via pest managers, the use of bed bug monitors or traps, quarantining of infested rooms and items, vacuuming, hot washing and drying infested items, the use of insecticides for small infestations, and the use of silicate products which are slow acting but have a long shelf life and a low possibility of resistance [7,8]. Examples of the advanced method include contained or circulated heat which might be integrated with some advanced chemical treatments. The use of dry heat and cold for bed bug control are relatively new and developing technologies, and the utility and effectiveness of these techniques has been subject to debate among pest experts [7,8,26]. We are interested in finding the optimal resource allocation that ensures the eradication of bed bugs in the system. In particular, we consider two scenarios.
In scenario I (S-I), the total number of exterminators, k i ,in each patch is fixed and we need to determine the number of exterminators that are using the advanced method in order to achieve the infestation free state (i.e., find optimal vector of (x 1 , x 2 , . . . , x n ) resulting in the minimum total cost subject to R 0 ≤ 1 ). Thus, the associated optimization problem reads as where R 0 is given by (5) with β ij and γ i given in (6) and (7), respectively, for i, j = 1, 2, . . . , n.
In scenario II (S-II), the total number of exterminators (K = n i=1 k i ) in all patches is fixed and we need to determine the optimal distribution of the exterminators to each patch and also the number of exterminators using the advanced method for reaching the infestation free state. Thus the optimization problem is 4. Examples. To give a specific example, we consider two cities in New Brunswick, Canada, namely, the city of Fredericton and the city of Saint John. Saint John is the most industrial city in New Brunswick and Fredericton is the capital city of the province. The distance between these two cities is about 113 km, and many people commute between two cities on a daily basis. Both cities are geographically close and there is not a significant difference between disinfestation standards and methods applied for the management of bed bug infestations. Given that Saint John has a history with bed bug outbreaks [4,5], the extensive human activities between these two cities put Fredericton at the risk of bed bug infestation as well. Thus, it is interesting to explore how disinfestation of bed bugs can be managed between these two cities. We use the 2−patch model as a prototype to explore the optimal resource allocation for disinfestation of bed bugs in this case. The basic reproduction number, R 0 , in equation (5) for two patches can be explicitly expressed as It is straightforward to note the following relation where R (i) 0 = β ii /γ i is the patch-specific basic reproduction number of patch i (i = 1, 2). This indicates that if R (i) 0 > 1 for i = 1 or i = 2, that is, if one patch undergoes an outbreak, then the bed bug infestation will spread throughout both patches.
The total cost (8) takes the following form where k 1 and k 2 are the number of exterminators in patches 1 and 2, respectively, x 1 and x 2 are the number of exterminators applying the advanced method in patches 1 and 2, respectively. For simplicity, we can assume that the time taken for an extermination method does not vary from one patch to another. That is, T 1 = T 11 = T 12 and T 2 = T 21 = T 22 , then the extermination rates for patches 1 and 2 are given by respectively. For simulation purpose, we fix the following parameters: T 1 = 60 , T 2 = 2 (days); p 11 = 0.2, p 12 = 0.7, p 21 = 0.4, p 22 = 0.6; c 11 = 4, c 12 = 2, c 21 = 3, c 22 = 4 (#contacts unit −1 day −1 ); and C 11 = 400, C 12 = 500, C 21 = 950, C 22 = 1000 (dollars).
For (S-I), we take k 1 = 8 and k 2 = 12. Thus, (9) simplifies to Minimize 0≤x1≤8,0≤x2≤12 For (S-II), we set K = k 1 + k 2 = 20 and hence (10) simplifies to Minimize 0≤x1≤k1≤20,0≤x2≤k2≤20 In optimization theory, the Lagrange multipliers method can be applied to find the local maximum or minimum of a continuous function subject to some equality constraints. For (S-I), we use the computer Algebra system Maple 14.00 (Maple 2010) to numerically minimize the total cost function in (15) subject to R 0 = 1, with two unknown variables x 1 and x 2 . The LagrangeM ultipliers() command in "Student[MultivariateCalculus]" subpackage, is then applied for this purpose. Note that according to our Theorem 2.1, the infestation diminishes to zero when R 0 ≤ 1, whereas the Lagrange multipliers method provides a real-valued pair (x 1 ,x 2 ) as the intersection point of the curves defined by the total cost function and R 0 = 1 in the x 1 − x 2 plane. This theoretical real-valued optimized pair is represented by the symbol • in Figure 2. A simple search procedure is then implemented to locate the practical optimal pair (x * 1 , x * 2 ) among the neighboring feasible integer pairs to obtain the minimum cost. Note that the cost function (13) is an increasing function of x 1 and x 2 , thus, increasing the number of exterminators applying the advanced method results in an increase of the total cost. This property of the cost function implies that in the x 1 − x 2 plane, only three integer pairs: ( x 1 , x 2 ), ( x 1 , x 2 ), and ( x 1 , x 2 ) need to be checked (See Figure 2). Here for a real-valued number x, x denotes the largest integer no larger than x and x is defined to be the smallest integer no smaller than x. Thus, the practical optimal pair will be the one with lowest total cost value. If the cost was identical for two pairs, the one with the smallest reproduction number R 0 would be chosen. The latter is due to the fact that a system with smaller R 0 takes less time to be disinfested.
For example, with the parameters chosen as above, the theoretical optimized pair of the problem in (15) is (x 1 ,x 2 ) = (3.9, 7.4) represented by symbol • in Figure  2. The three pairs to be checked are: (3,8), (4,7), and (4,8). For the first two, R 0 > 1, thus the practical optimal pair is (x * 1 , x * 2 ) = (4, 8) with R 0 ≈ 0.95 and the total cost $15, 400. Theoretical and practical optimal solutions and a non-optimal feasible solution to (S-I). The dashed line represents the cost function in (13) with k 1 = 8 and k 2 = 12. The solid curve is defined by R 0 = 1. The symbol * indicates the practical optimal pair (4,8) and the circle represents a non-optimal feasible pair (6, 10).
To find the optimal solution for (S-II), the procedure employed in (S-I) is looped for k 1 from 1 to k 1 +k 2 = K. The optimal solution is then given by the combination (k 1 , k 2 , x 1 , x 2 )that gives the minimum cost among all feasible combinations. For the above chosen parameters, the optimal combination is (k * 1 , k * 2 , x * 1 , x * 2 ) = (12,8,4,8). This indicates that to achieve the infestation free state with the minimum cost, among 20 available exterminators, 12 should be allocated to patch 1, of which 4 should apply the advanced method, and the rest 8 exterminators should be allocated to patch 2, all of them should apply the advanced method. The optimal cost is $15, 000 and the corresponding R 0 ≈ 0.94.
We define a proportion-index, I e , to be the indicator of successful management. That is, if I i < I e for i = 1, 2, . . . , n, then we say the infestation can be negligible. We use t e to denote the index-time at which the indicator state is achieved. In our simulations, we take I e = 0.01 (i.e, if the infested units in a patch are less than 1%, then the indicator state is achieved).
Based on the 2011 Canada census profile provided in [29], the total number of private dwellings in Fredericton and Saint John were 41, 581 and 84, 542, respectively, which can be regarded as reasonable estimates of the total number of units in these two cities, namely, N 1 and N 2 , in system (1).
Suppose, initially, there were 2, 079 infested units in Fredericton, and 8, 454 infested units in Saint John. This gives the initial value of system (4): I 1 (0) = 0.05 and I 2 (0) = 0.1. Assume that there are 8 exterminators in Fredericton and 12 exterminators available in Saint John, i.e.,k 1 = 8, k 2 = 12. This corresponds to scenario (S-I) and it follows from our proposed optimization procedure that the optimal pair is (x * 1 , x * 2 ) = (4,8), and the corresponding basic reproduction number is R 0 ≈ 0.94, while the minimum total cost is C = $15, 400, and the index-time taken to reduce the infested units in both patches to be below 1% is t e = 9.2 days; while at a non-optimal pair, say, (x 1 , x 2 ) = (6, 10), the corresponding R 0 ≈ 0.76 and the total cost is C = $17, 500. It takes t e = 1.7 days to reach the indicator state, but at the cost of extra $2, 100 (See Figure 3). In practice, there should be a trade-off between management time and the corresponding cost and more information should be collected for the decision makers to determine the best strategy.
The procedure of finding the optimal pair that gives practical resource allocation strategy we proposed for the case with n = 2 can be extended to general case with n > 2 under similar management scenarios as (S-I) and (S-II). However, in practice, it is very difficult to express R 0 defined in (5) explicitly as a function of transmission and extermination rates when n > 2. In this case, we employ the idea of the method of exhaustion and loop n times through x i s from 1 to k i to numerically compute the corresponding R 0 . This will allow us to obtain the cost associated to each n-dimensional vector (x 1 , x 2 , . . . , x n ) satisfying R 0 < 1 . Finally, a simple search among all feasible solution vectors will find the minimum cost. Note that the minimum cost can be achieved at multiple feasible solution vectors. The optimal solution vector is the one that gives the minimum cost with the smallest corresponding R 0 . For example, for (S-I) with n = 4 and k 1 = 8, k 2 = 12, k 3 = 12, k 4 = 20, Assume that the extermination time taken for the preliminary and the advanced methods are T 1 = 60 and T 2 = 2 (days), respectively and the associated costs for both methods in those 4 patches are C 11 = 400, C 12 = 500, C 13 = 450, C 14 = 550, C 21 = 950, C 22 = 1000, C 23 = 1000, C 24 = 1100 (dollars). Then the optimal solution vector is (x * 1 , x * 2 , x * 3 , x * 4 ) = (8,12,11,9), which gives the optimal cost $47, 000 with R 0 ≈ 0.9947. 5. Summary and discussion. In this paper we have proposed a modeling framework for studying the spread of bed bug infestation among multiple patches. Susceptible units in each patch can become infested due to human activities among connected patches. Assuming irreducible transmission matrix F = (β ij ) and constant extermination rates, we have shown that our model exhibits global threshold dynamics in the sense that if the basic reproduction number is less than 1, then the infestation free equilibrium is globally asymptotically stable and all involved patches will eventually become bed bug infestation free. Further, if the basic reproduction number is larger than 1, then the positive equilibrium is globally asymptotically stable and the bed bug infestation persists in all involved patches.
To manage the bed bug infestation, by considering costs associated with extermination processes, we set up the corresponding optimization problems. For the two practical scenarios considered here, we proposed procedures to find the optimal resource allocation strategies. The optimal resource allocation for disinfestation management we explored is optimal in the sense that the total cost is minimum subject to R 0 ≤ 1. This does not imply that the associated R 0 is the smallest among all feasible resource allocations. Indeed, R 0 obtained at any other feasible resource allocations is always smaller, resulting in a quicker disinfestation process, but a higher cost.
Note that in our modelling framework, it is assumed that there is no time delay for a susceptible unit to become infested after encountering the infested units via human activities. It is also possible to incorporate the time delay into the modelling and follow the approaches employed in [24,30] to establish the global threshold dynamics. In addition, the idea can be extended to the case with non-constant transmission rates and Allee effect [17], to the case with varying population sizes [20] and to the case with nonlinear incidence rate [19,36].
For our model, let E = n i=1 [0, 1], then it is easy to verify that E is positively invariant with respect to system (4). Note that for i = j and I ∈ E, (Df ) ij (I) = ∂f i (I)/∂I j = β ij (1 − I i ) ≥ 0.
Thus the map f defined by the righthand side of (4) is continuously differentiable and is cooperative and Df (I) is irreducible. That is, condition (1) of Lemma A.1 is satisfied. Since I 0 = 0 is the IFE, and E is positively invariant, condition (2) of Lemma A.1 is easily met. For any λ ∈ (0, 1) and I ∈ E with I i > 0 for i = 1, . . . , n, we have f is strictly sublinear and hence condition (3) of Lemma A.1 is verified. Consequently, we arrive at the threshold dynamics stated in Theorem 2.1.