DYNAMICS OF A METAPOPULATION EPIDEMIC MODEL WITH LOCALIZED CULLING

. A two–patches metapopulation mathematical model, describing the dynamics of Susceptibles and Infected in wildlife diseases, is presented. The two patches are identical in absence of control, and culling activities are performed in only one of them. Firstly, the dynamics of the system in ab- sence of control is investigated. Then, two types of localized culling strategies (proactive end reactive) are considered. The proactive control is modeled by a constant culling eﬀort, and for the ensuing model the disease free equilibrium is characterized and existence of the endemic equilibrium is discussed in terms of a suitable control reproduction number. The localized reactive control is modeled by a piecewise constant culling eﬀort function, that intro- duces an extra–mortality when the number of infected individuals in the patch overcomes a given threshold. The reactive control is then analytically and numerically investigated in the frame of Filippov systems . We ﬁnd that localized culling may be ineﬀective in controlling diseases in wild populations when the infection aﬀects host fecundity in addition to host mortality, even leading to unexpected increases in the number of infected individuals in the nearby areas.

(Communicated by the associate editor name) Abstract. A two-patches metapopulation mathematical model, describing the dynamics of Susceptibles and Infected in wildlife diseases, is presented. The two patches are identical in absence of control, and culling activities are performed in only one of them. Firstly, the dynamics of the system in absence of control is investigated. Then, two types of localized culling strategies (proactive end reactive) are considered. The proactive control is modeled by a constant culling effort, and for the ensuing model the disease free equilibrium is characterized and existence of the endemic equilibrium is discussed in terms of a suitable control reproduction number. The localized reactive control is modeled by a piecewise constant culling effort function, that introduces an extra-mortality when the number of infected individuals in the patch overcomes a given threshold. The reactive control is then analytically and numerically investigated in the frame of Filippov systems.
We find that localized culling may be ineffective in controlling diseases in wild populations when the infection affects host fecundity in addition to host mortality, even leading to unexpected increases in the number of infected individuals in the nearby areas.
1. Introduction. In the last decades, emerging and re-emerging infectious diseases (ERIDs) have been responsible for significant economic and social impacts in both developed [25] and developing countries [34]. Jones et al. [33] showed that zoonotic infections originated in wildlife have accounted for the majority of ERID events since the 1940s and are representing a growing threat to global health. As a consequence, actions aiming at controlling diseases in wild populations can provide benefits on wildlife protection and conservation as well as on the sustainability of agriculture and public health.
However, disease control through vaccination and drug treatments, which represent common intervention measures in human and livestock infections, is often not feasible, or practical, in wildlife because of the lack of resources or the unavailability of suitable diagnostic tools [18]. In these circumstances, non-selective culling, which consists in the slaughtering of both infected and healthy individuals, represents the only available disease control strategy.
The effectiveness of culling relies on the assumption that, under a given threshold of host density, the population becomes too sparse, then, the number of potentially infectious contacts between infected and susceptible individuals become too low to allow the disease to spread and persist in the population [2]. Yet there are empirical evidences suggesting that culling has been ineffective in reducing disease burden in wildlife populations for different infections, such as rabies in canids [37,46], facial tumor disease in Tasmanian devil [36], and bovine tuberculosis in European badger [27].
The main causes for the failure of culling in the eradication of diseases in wild populations have been ascribed to compensatory mechanisms, such as the densitydependent positive feedbacks on recruitment and dispersal triggered by the excess of mortality due to culling [12,43]. Specifically, by reducing the host population abundance, culling can reduce the density-dependent constrains on host birth rate, thereby producing a flush in new susceptible individuals in the population [29,46]. These new susceptibles represent a reservoir for new infections, which nullifies the expected benefits of disease control campaigns or, in some cases, even increases the disease burden in the population [5,7,15] or the duration of the epidemic [9,10].
On the other hand, different studies showed that culling may disrupt the host social structure increasing the animal home range and prompting long-distance movement and dispersal (perturbation hypothesis), thus increasing the probability of potentially infectious contacts between neighboring groups [13,39,42]. In particular, studies carried in the British Isles on the effectiveness of localized proactive and reactive culling as measures for bovine tuberculosis (Mycobacterium bovis) control in European badger showed that culling led to a decrease in disease burden in the controlled lands, while caused an increase in the infections in the nearby areas [3,4,26,45]. As expected from the perturbation hypothesis, the increase of M. bovis prevalence was associated with expanded home ranges and more frequent migration events in badger [4].
In this work, we will consider a metapopulation epidemic mathematical model, described by a system of ODEs for Susceptibles and Infected, with two patches which are identical in absence of control and localized culling control which is performed in only one of them. We will show by means of qualitative analysis and numerical simulations that, also in the absence of the density-dependent compensatory effects on recruitment and dispersal, localized culling may be ineffective in controlling diseases in wild populations and it may even lead to unexpected increases in the number of infected individuals in the nearby areas.
Here, we investigate the effectiveness of localized culling using a metapopulation Susceptibles-Infected (SI) model with density-dependent mortality of the host and disease-induced host sterilization. Disease-induced fecundity reduction (or sterilization) of the host has been frequently observed in host-pathogen interactions. Such kind of disease-induced fecundity reduction of the host has been frequently observed in host-pathogen interactions. Mathematical models describing the dynamics of fecundity-reducing infections have been developed for rabies -where the fecundity reduction is generated by the severely debilitating nature of the infection and the frenzied behavior induced -in different mammal species [6], cowpox in wild rodents [40], and several infections of invertebrate hosts, such as crustaceans and mollusks [11]. Moreover, we will show that the counter-intuitive effects of culling hold also relaxing the hypotheses of disease-induced host sterility and densitydependent host mortality in the infection model.
2. The metapopulation mathematical model. We expand a traditional SI epidemic framework developed for host-sterilizing infections in wild animals to describe the effect of localized proactive and reactive culling on the control of wildlife diseases in a two-patches metapopulation model (see e.g. [2,5,19]). From an epidemiological point of view, hosts can be subdivided into compartments with respect to the infection: the susceptible compartments S j (t) -i.e., healthy individuals that can be infected by the pathogen -and the infected compartments I j (t) -i.e., diseased individuals that can infect other individuals; subscript j(= 1, 2) defines the patch in the metapopulation to which the epidemiological compartment belongs.
In order to describe localized culling control in the model we assume that culling activities can be performed in one of the two patches only (specifically patch 1), while the other one (patch 2) is uncontrolled. We also assume the two patches are identical in the absence of control, namely the demographic and epidemiological parameters are the same.
Then, the metapopulation model can be represented by the following set of four ordinary differential equationṡ where the upper dot is used to denote the time derivative. N j = S j + I j represents the total host population in patch j. All the parameters are positive constants: r = ν − µ represents the intrinsic growth rate, being ν the fertility rate and µ the natural mortality rate; γ and α represent the density-dependent mortality rate of a disease-free host population and the additional mortality rate caused by the infection, respectively; and β represents the transmission rate between infected and susceptible individuals. Moreover, we define D to be the per-capita dispersal rate from a patch to the other; and c to be the culling effort. In order to analyze the effects of both localized proactive and reactive culling strategies in model (1), we implement two different control functions. Thus, the localized proactive culling strategy in model (1) is obtained by imposing culling effort c to be a constant function and the localized reactive culling strategy by imposing culling effort c to be a piecewise constant function as in [20] c =c 2 wherec represents the maximum effort enforceable, and θ represents the threshold for the detection of the infection. Expression (3) implies that population N 1 in patch 1 undergoes an extra-mortalityc due to culling when the number of infected individuals in the patch is higher than threshold θ, while no culling activities are performed when the infection is under the threshold of detection. Therefore, threshold θ introduces a discontinuity, so that model (1) with (3) becomes a discontinuous piecewise-smooth system (also called Filippov system) in which sliding motions are possible on the manifold separating the region (in state space) where control is allowed from that where it is not [28]. More precisely, by introducing the so-called switching manifold and denoting with f (1) (resp. f (2) ) the RHS of model (1) with c = 0 (resp. c =c), then sliding occurs on the sliding set namely, where the components of two vector fields f (i) , i = 1, 2, transversal to Σ are 'pushing' in opposite directions, forcing the state of the system to remain on the switching manifold and slide on it. Σ s terminates in Σ when a tangency occurs, i.e. when f vanish, implying that f (1) or f (2) are tangent to the switching manifold. Tangencies are strategically important for bifurcation analysis. As first pointed out by Filippov, sliding motions obey the ODE system having as RHS the unique convex combination of f (1) and f (2) parallel to Σ s [28], i.e., Equilibrium points of system (5) are called pseudo-equilibria for model (1)-(3) and correspond to a stationary sliding solution. The analysis of planar piecewise-smooth autonomous systems has been widely developed in recent years, see e.g. [31,35]. However, there is still no obvious classification of the dynamics of higher dimensional piecewise-smooth systems [30]. Then, the behaviors exhibited by model (1) with localized reactive culling (3), which represents a four-dimension piecewise-smooth system, are too complex to be fully investigated analytically. As a consequence, we numerically investigate system (1)-(3) through simulations performed with the event-driven method developed by Piiroinen and Kuznetsov [38] for sliding systems to find model attractors. A bifurcation analysis in the parameter space identified by the reactive control parameters,c and θ, has been also performed with SLIDECONT [21], which is a software based on the package AUTO [24] to continue solutions to nonlinear boundary-value problems via orthogonal collocation.
3. Results. It can be shown that model (1) is consistent, namely its state variables remain positive for any trajectory starting from positive conditions, as stated by the following theorem: Theorem 3.1. If we consider positive initial conditions S 1 (0), I 1 (0), S 2 (0) and I 2 (0), then the solutions of the differential system (1) are positive at each time t > 0.

DYNAMICS OF A METAPOPULATION EPIDEMIC MODEL WITH LOCALIZED CULLING 5
Proof. We shall prove the statement by contradiction. Let us start by considering the state variables S 1 (t) and S 2 (t) and let t 1 > 0 be the first time instant when S 1 (t)S 2 (t) = 0. Since the initial conditions are positive, the variables S 1 and S 2 (and hence their product) are positive in [0, t 1 ).
Thus, the existence of t 1 is definitely excluded.
With the same arguments one can prove the positivity of the infected state variables I 1 (t) and I 2 (t).
3.1. Metapopulation model without control. Before exploring the effects of localized proactive and reactive culling on the infection dynamics, we analyze the main features of model (1), deriving the equilibria in the absence of control (c = 0) and the expression for the basic reproduction number, R 0 (which corresponds to the average number of secondary infections caused by a single primary infection in a totally susceptible population at the disease-free equilibrium). In the absence of control (c = 0), model (1) assumes a symmetric form, which allows us to prove that: (1), then at any equilibrium point the size of the susceptible (resp. infected) compartment in patch 1 is equal to the size of the susceptible (resp. infected) compartment in patch 2.
Proof. Let us consider the algebraic system obtained by setting the RHS of model (1), with c = 0, equal to zero, namely: By contradiction, we assume that there exists a positive solution to system (6), say Then, two cases must be considered: which reduce toN 1 <N 2 and a contradiction arises. A similar argument applies ifS 2 >S 1 . •Ī 1 =Ī 2 In order to fix the ideas, let us assumeĪ 1 >Ī 2 (the opposite case follows in analogous way). IfS 1 >S 2 , then, analogously to what happens ifĪ 1 =Ī 2 , the inequalities (7) reduce to γN 1 +βĪ 1 < γN 2 +βĪ 2 , which is in contrast with the assumptions Otherwise, if 0 <S 1 ≤S 2 , then, the inequalities (7) must be reversed, namely yielding: On the other hand, by handling the (6b) and (6d), one obtains which, for the positivity of theĪ i , i = 1, 2, reduce to BeingS 1 ≤S 2 , the last inequality implies thatN 1 <N 2 . However, adding the (9) and (11) leads to a contradiction:N 1 >N 2 .
In conclusion, ifS 1 =S 2 = 0, then the (10) are still valid and reduce tō I 1 <Ī 2 , which is again in contrast with the initial assumption.
In other words, in searching model equilibria, it is not restrictive to assume the two patches are identical. In particular, model (1) with c = 0 has always a trivial equilibrium, E 0 = [0, 0, 0, 0], and a disease-free equilibrium, represents the carrying capacity for a disease-free host population. It is easy to prove that E 0 is always unstable if r > 0, while E 1 is asymptotically stable if and only if the basic reproduction number of model (1), R 0 , is lower than 1 [1]. The basic reproduction number can be calculated as the spectral radius of the next generation matrix, R 0 = ρ(F V −1 ), where F and V are defined as Jacobian matrices of the new infections appearance and the other rates of transfer, respectively, calculated for infected compartments at model (1) disease-free equilibrium [22,44]. Then, model (1) basic reproduction number can be defined as with K given in (12) and coincides with that for the corresponding homogeneous mixing model. When R 0 > 1 in equation (13), the disease-free equilibrium is unstable and there exists a unique asymptotically stable positive endemic equilibrium, We notice that expressions (14), defining the endemic equilibrium in (1) for c = 0, correspond to the endemic equilibrium for the associated homogeneous mixing model (see [2], under the assumption of negligible incubation period) and stability properties follow straightforwardly. Note also that necessary condition for beingĪ positive is that hence we assume it always fulfilled in the following.

DYNAMICS OF A METAPOPULATION EPIDEMIC MODEL WITH LOCALIZED CULLING 7
3.2. The effect of localized proactive culling. Analyzing metapopulation model (1)-(2) -i.e., with localized proactive culling strategy -, we find that there always exists a trivial equilibrium, Ec 0 = [0, 0, 0, 0]. Conversely than the case with c = 0, the trivial equilibrium Ec 0 is not always unstable when r > 0. Indeed, linearization of system (1)-(2) around Ec 0 leads to Ec 0 locally asymptotically stable if and only if as can be proved by Sylvester criterion [32], being the Jacobian matrix of system (1)-(2) evaluated at Ec 0 a symmetric matrix. Expression (16) suggests that, for frequent dispersers (i.e., species with D > r, as defined in [14]), sufficiently high levels of constant culling efforts (c >ĉ) in one of the patches can lead to the extinction of the entire population.
When expression (16) is not satisfied, the trivial equilibrium is unstable and there exists a positive disease-free equilibrium, Ec 1 , as stated by the following theorem: (16) is not verified, then model (1)-(2) admits an unique positive disease-free equilibrium Proof. Denote the generic disease-free equilibrium equilibrium of model (1)-(2) by The componentsK 1 ,K 2 are the solutions of the algebraic system obtained by setting the RHS of equations (1a) and (1c), with c =c, equal to zero, namely: (r − D)K 2 − γK 2 2 + DK 1 = 0, yielding the admissibility condition: Substituting the expression ofK 1 into (17), one obtainsK 2 = 0 as solution of with , and Let us differentiate the equation (19), yielding: If ∆ ≤ 0, then f is an increasing function. Otherwise, if ∆ > 0, then (3.3) admits two real solutions: that are relative minimum/maximum points for f and it is straightforward to check thatK In any case, we have Thus, exactly one root of (19) satisfies the admissibility condition (18). If D > r, it is not obvious that such a root is positive: one can easily check that (16) is (resp. is not) fulfilled. Hence, in summary, system (1)-(2) admits an unique disease-free equilibrium Ec 1 if condition (16) is not satisfied; there are none if (16) is verified.
Thus, in the presence of localized proactive culling the disease-free equilibrium becomes Ec 1 = [K 1 , 0,K 2 , 0], whereK 1 andK 2 can not be easily written explicitly, since they come from solutions of the third-order equation (19).
When proactive culling as in (2) is implemented in model (1), we can compute through the next generation matrix method the control reproduction number, R C , similarly to R 0 . The control reproduction number is defined as the average number of secondary infections produced by a single infected individual in a susceptible population at its disease-free equilibrium experiencing culling effortc. Its expression is with (for details see Appendix A). When it exists, disease-free equilibrium Ec 1 is asymptotically stable if and only if the control reproduction number R C in (20) is lower than 1; on the other hand, if R C > 1 the disease-free equilibrium is unstable [22,44]. As far as the existence and the number of model endemic equilibria, an analytical investigation is very hard to perform. Therefore, numerical analyses are carried out to understand the effect of localized proactive culling on the disease dynamics. Extensive numerical tests suggest that there are no endemic equilibria when R C < 1; otherwise, when R C > 1, there exists an unique asymptotically stable endemic equilibrium, say Ec 2 = [Ŝ 1 ,Î 1 ,Ŝ 2 ,Î 2 ].
In Fig. 1 we show the effect of the host dispersal rate (D) and the localized proactive culling effort (c) on the existence and stability of model (1)-(2) equilibria through bifurcation analysis [41]. The curve T C 0 defines a transcritical bifurcation representing the threshold for host extinction (in (16)) that separates the region in which the total host population goes extinct from the region in which model (1)-(2) converges toward a disease-free equilibrium. The curve T C 1 defines a transcritical bifurcation representing the threshold for infection establishment (R C = 1) that separates the region in which the pathogen fails to establish itself (disease-free equilibrium) from the region in which the pathogen is able to invade the host population and model (1)-(2) converges toward an endemic equilibrium.  The curves T C 0 (i.e.,c =ĉ in (16)) and T C 1 (i.e., R C = 1), which represent transcritical bifurcations, delimit three different regions in the parameter space [D,c] where system (1), in the presence of localized proactive culling (2), converges to host extinction, disease-free equilibrium, or endemic equilibrium. Parameter values for model (1)- (2) have been fixed to: r = 0.9; µ = 0.2; K = 600; α = 0; R 0 = 10 (or rather β = 0.01833, see (13)).
Proof. It is an immediate consequence of equations (24a) and (24b).
From expression (23), we notice that necessary conditions for localized proactive culling to increase (instead of decrease) the number of infected individuals in patch 2 are: 2D < r (i.e., hosts are very infrequent dispersers) and R 0 > 2 (i.e., in the corresponding uncontrolled model the first infected individual can, on average, infect more than two individuals).

Remark 1.
In the case the total host population at endemic equilibrium for c = 0 (15) is much lower than the carrying capacity K, the effect of the density-dependent mortality rate at the endemic equilibrium is negligible, which corresponds to assume γ → 0. In this scenario, it is straightforward to check that condition (23) for beinĝ I 2 (0) > 0 reduces to: 2D < r. Instead, if, in addition to setting γ → 0, we also relax the hypothesis of disease-induced sterility (see Appendix B), then such a condition becomes slightly more complex: where α > r (which represents the necessary condition for the existence of the endemic equilibrium in the absence of control, see Appendix B). However, once again, formula (30) indicates that having 2D < r is necessary for proactive culling to increase the endemic value of infected in the uncontrolled patch.
To investigate more deeply the counter-intuitive result obtained in Theorem 3.6, we numerically compute the parameter conditions for which the number of infected individuals in patch 2 at the endemic equilibrium is higher in the presence of localized proactive control than with the do-nothing alternative (i.e.,Î 2 (c) >Ī). The numerical analyses presented here are performed by exploring the effects of different values of host dispersal (D), proactive culling effort (c), and disease basic reproduction number of model (1) with c = 0 (R 0 , given by (13)) on control effectiveness and by keeping the host demographic parameters constant in the simulations. Specifically, we set the host demographic parameters as in [5].
In Fig. 2 we show the parametric regions where conditionÎ 2 (c) >Ī is satisfied for different values of basic reproduction number of the corresponding uncontrolled model (R 0 ) in the parameter space [D,c]. Fig. 2 highlights that for low values of host dispersal rates (D), localized culling is ineffective in reducing the infection outside the control zone (i.e.,Î 2 (c) >Ī) for a wide range of culling rate values, also when the basic reproduction number is relatively low.
The dynamics of infected individuals in patch 2 corresponding to the parameter space along the thin dotted line in Fig. 2 is illustrated in Fig. 3. In particular, in  Fig. 3. Unspecified parameters as in Fig. 1. 3.3. The effect of localized reactive culling. Let us consider model (1) with the piecewise constant function (3) representing the culling effort in the localized reactive control strategy. We are interested in investigating the long term dynamics of model solutions and, in particular, of sliding motions, occurring on the sliding set (4). To this aim, the following results concerning the existence and stability of stationary sliding solutions (also called pseudo-equilibria) can be proved: Theorem 3.7. Necessary condition for the existence of pseudo-equilibria for model whereÎ 1 (resp.Ī) is the number of infected individuals in patch 1 at endemic equilibrium Ec 2 (resp. E 2 ) of model (1) with c =c (resp. c = 0).
Since in any case we find a contradiction, θ must belong to the interval (31).
Proof. From (32), the Jacobian matrix of system (5) at a pseudo-equilibrium reads with d 1 , d 2 and d 3 given in (37). With simple algebraic calculations, we derive the characteristic polynomial of J (say, P (λ)): Being l i > 0, ∀i = 0, . . . , 2, the presence of positive real roots for P (λ) is excluded in virtue of Descartes' rule of sign. According to Routh-Hurwitz theorem, also complex roots with positive real part are not admissible if, and only if, l 1 l 2 − l 0 > 0, which corresponds to (36).
The analytical condition (36) provided in Theorem 3.8 depends in a complicated manner on unknown pseudo-equilibrium values as well as on crucial parameters, like β, γ, D, θ; hence, it is difficult to give easier sufficient conditions for its positivity. However, in the specific case the hosts are very infrequent dispersers (i.e. D → 0), which corresponds to the scenario where culling leads to an increase of infected individuals in the uncontrolled patch, expression (36) is always verified. In this case, a pseudo-equilibrium of model (1)-(3) is locally asymptotically stable.
Since the analytical results provided in Theorems 3.7 and 3.8 are not exhaustive, we perform different sets of numerical simulations on a wide range of parameters combinations for θ andc (as in Fig. 4) and for different initial conditions of model (1) variables. For each pair of parameters θ andc, we find a unique attractor regardless of the initial conditions chosen for the simulations. Specifically, we find that the only suitable attractors for model (1)-(3) are either stable equilibria or pseudo-equilibria (which represent points of the sliding set (4).
The bifurcation analysis of epidemic model (1)-(3) in the parameter space [θ,c], derived from the continuation of stable equilibria and pseudo-equilibria found in the numerical simulations, is illustrated in Fig. 4. Fig. 4 shows that, for large values of culling threshold (θ), model (1) equilibrium (independent from the values assumed byc) corresponding to endemic equilibrium (14) as defined in continuous model (1) in the absence of control, namely E 2 (see region 1 in Fig. 4). By decreasing threshold θ, equilibrium E 2 found in region 1 undergoes a boundary-node bifurcation for θ = I 0,0 Fig. 4), where I θ,c j denotes the steady-state value of infected individuals in patches j for localized reactive control (3) with parameters θ,c. By crossing BN1 and entering region 2 the original E 2 equilibrium disappears and a pseudoequilibrium characterized by I θ,c 1 = θ appears (see region 2 in Fig. 4). Curve BN2 represents a boundary-node bifurcation: by crossing it and entering region 3, the pseudo-equilibrium characterized by I θ,c 1 = θ disappears and a stable equilibrium characterized by I θ,c 1 > θ appears. Boundary-node bifurcations BN1 and BN2 correspond to the vanishing of the vector fields f  Fig. 6). 4. Discussion. In this paper, we analyze the effect of localized proactive and reactive culling on the disease dynamics in a metapopulation model with two patches (one with control and the other one without control). Proactive culling is described through a classical ODE continuous system, while reactive culling is described through a discontinuous piecewise-smooth system where the control activities are implemented when the number of infected individuals exceeds a given threshold. Models implementing sliding control have been already developed in recent years for different infections, such as West Nile Virus [48], avian influenza [17,16], and SARS [47]. Here, we find that, localized culling implemented in one of the patches  2), we provide the necessary conditions for this counter-intuitive outcome to occur (see Theorem 3.6). In details, the equilibrium value of infected individualsÎ 2 in presence of proactive culling can increase only if 2D < r (hosts are very infrequent dispersers) and R 0 > 2, namely if the basic reproduction number of the corresponding uncontrolled model is at least twice the endemicity threshold. In addition, we numerically find that the number of infected individuals in the uncontrolled patch eventually peaks for intermediate values of culling and then decreases only for high level of culling effort (see Fig. 3). On the other hand, in sliding model (1)-(3), we numerically find that, for intermediate levels of disease detection, localized reactive culling increases the infection burden in the uncontrolled patch regardless of the culling effort applied in the disease control (see Figs. 4 and 6B). The biological explanation for this unexpected effect relies on the remark that, when the dispersal rate (D) is low, the introduction of culling induces an increase in the number of susceptible individuals in both the controlled and uncontrolled patches (see Theorem 3.6 and Proposition 1). This leads to a flush of new susceptibles entering the uncontrolled patch, then providing a bust to the infection transmission. We have proved the generality of this result by showing that it does not depend on the specific assumptions made in model (1), but it holds also relaxing the hypotheses of disease-induced host sterility and density-dependent mortality in the infection model (see Remark 1). The epidemiological and economic conditions leading to the ineffectiveness of disease control through culling have also been investigated in optimal control frameworks. Bolzoni et al. [8] showed that reactive culling implemented around the peak of infection represents an optimal control strategy only when the basic reproduction number of the infection is low and the costs of control are high.
In conclusion, this paper shows that also in the absence of the density-dependent compensatory effects -which were previously associated to disease control failure -, localized culling may represent an ineffective strategy in limiting infectious diseases in wildlife.
Appendix A. The control reproduction number R C . Following the procedure and the notations in [22,44], we prove that the control reproduction number of model (1)-(2), R C , is given by (20).