Model-Based Evaluation of Strategies to Control Brucellosis in China

Brucellosis, the most common zoonotic disease worldwide, represents a great threat to animal husbandry with the potential to cause enormous economic losses. Brucellosis has become a major public health problem in China, and the number of human brucellosis cases has increased dramatically in recent years. In order to evaluate different intervention strategies to curb brucellosis transmission in China, a novel mathematical model with a general indirect transmission incidence rate was presented. By comparing the results of three models using national human disease data and 11 provinces with high case numbers, the best fitted model with standard incidence was used to investigate the potential for future outbreaks. Estimated basic reproduction numbers were highly heterogeneous, varying widely among provinces. The local basic reproduction numbers of provinces with an obvious increase in incidence were much larger than the average for the country as a whole, suggesting that environment-to-individual transmission was more common than individual-to-individual transmission. We concluded that brucellosis can be controlled through increasing animal vaccination rates, environment disinfection frequency, or elimination rates of infected animals. Our finding suggests that a combination of animal vaccination, environment disinfection, and elimination of infected animals will be necessary to ensure cost-effective control for brucellosis.


Model parameters and values
Constant parameters. Table S1 shows the epidemiological parameter values for brucellosis in literature. The unit of the model parameters is one year. The incubation period of human brucellosis is about two weeks [1], so the clinical outcome rate for exposed people is σ = 26 per year. From the China Statistical Yearbook [2] and the China Animal Husbandry Statistical Yearbook [3], one can obtain the demographic parameter values (including human populations birth rate b h and death rates d h , sheep recruitment and slaughter rate b) for mainland China and the 11 provinces with the highest incidence for human brucellosis listed in Table S2.  Estimated values of transmission rates. We fix the human indirect transmission rate β hw , and assume that β hw = 0.5 for mainland China and 11 selected provinces with the highest incidence for human brucellosis. By using least-squares fitting method in DEDiscover software, the estimates of β s , β sw and β h for these 11 provinces and mainland China can be obtained, and are shown in Table S3.

Dynamical behavior
For the brucellosis model: As the last four equations are independent of the first three equations for system (1), we only need to consider the following system: For the general incidence rate β sw S s f (W ), nonnegative functions f (W ) is assumed to be differentiable with W , and thus solutions to system (2) with nonnegative initial conditions exist and are unique. Throughout we also assume the following properties of function f (W ), which are biologically reasonable: W is monotone nonincreasing. Notice that mass action incidence, saturating incidence, and standard incidence of system (2) satisfy assumptions (H1)-(H3). Omega limit sets of system (2) are contained in the following bounded region in the non-negative cone of R 3 : Region X is positively invariant with respect to system (2). It is obvious that any solution of system (2) with nonnegative initial values is nonnegative and system (2) has one disease-free equilibrium P 0 = (N s , 0, 0). We derive the basic reproduction number of system (2) by the next generation matrix formulated in Diekmann et al. [5,6]. We order the infection variables first by disease state, only needing the vector x = (I s , W ) T . Considering the following auxiliary system: Following the recipe from van den Driessche and Watmough [7] to obtain: is the derivative of f (W ) with respect to W at disease-free equilibrium. The basic reproduction number is defined as the spectral radius of the nonnegative matrix F V −1 ; it is easy to obtain where R e 0 = kβswNsf (0) bδ and R i 0 = βs b are partial reproductive numbers due to environment-toindividual transmission and individual-to-individual transmission, respectively.
Stability of the disease-free equilibrium. System (2) is a cooperation system, and the following theorem shows that the disease-free equilibrium of system (2) is globally asymptotically stable in the region X.
Theorem 2.1. Supposing that assumptions (H1)-(H3) hold, then the following conclusions hold for the system (2). (a) If R 0 < 1, then the disease-free equilibrium P 0 of system (2) is globally asymptotically stable in the region X. (b) If R 0 > 1, then the disease-free equilibrium P 0 of system (2) is unstable.
Proof. One can obtain the following equation by using assumptions (H1) and (H3), which also means f (0)W ≥ f (W ). Hence, for system (3), it is easy to obtain: Define the Lyapunov function: Then the derivative of L along the system (3) is: If R 0 < 1, dL 1 dt = 0 implies that b T x = 0, thus I s = 0, W = 0. Therefore, the largest invariant set of Ψ is the singleton P 0 . By LaSalle's invariance principle [8], P 0 is globally asymptotically stable in the region X when R 0 < 1.
If R 0 > 1 and x > 0, it follows that: There must exist dL 1 dt > 0 in a small enough neighborhood of P 0 in the interior of X. Therefore, solutions in the interior of X sufficiently close to P 0 move away from P 0 provided R 0 > 1, and thus P 0 is unstable. The proof is end.
The existence and stability of the endemic equilibrium. Firstly, we will show the existence of the endemic equilibrium of system (2). It is assumed that P * = (S * s , I * s , W * ) is an endemic equilibrium of system (2), and satisfies the following equilibrium equations: We have W * = kI * s δ and S * s = N s − I * s . Therefore, the equilibrium of system (2) is equal to the following system: Equation (7) has a zero solution I * s = 0. Let, For equation (8), we can obtain: Due to F (0) = 0 and F (N s ) = −bN s < 0, so the sufficient and necessary condition for the existence of positive equilibrium of system (7) is: Hence, we can conclude that for system (2) there at least exists an endemic equilibrium if R 0 > 1.
Theorem 2.2. Suppose that assumptions (H1)-(H3) hold. If R 0 > 1, then the endemic equilibrium P * of system (2) is unique and globally asymptotically stable in the interior of X.
Proof. Let h(a) = 1 − a + ln a for all a > 0, then it is easy to verify that: Differentiating D 1 and D 2 with t along solution curves of system (2) and using equilibrium equations (6) to simplify, one can obtain: Through using inequality (9), we can obtain the second and third inequality with the following equations: Following from assumptions (H2) and (H3), we have Hence, one can obtain: Similarly, we have the following equation: Hence, we can define the following Lyapunov function: It follows that: Moreover, the equality dL 2 dt = 0 holds if and only if S s = S * s , I s = I * s and W = W * . Thus, P * is the only invariant set of system (2) in {(S s , I s , W ) ∈ X : dL 2 dt = 0}. By LaSalle's invariance principle [8], P * is globally asymptotically stable and unique in the interior of X when R 0 > 1. This completes the proof.
The endemic equilibrium of human population. From previous analysis we know that the endemic equilibrium of system (2) is unique and globally asymptotically stable if R 0 > 1.
(1) The indirect transmission rate of human population is standard incidence, and b h > d h . When system (2) is stable, the last four equations of system (1) will become the following equations: where We make the following scaling transforms for system (10), and we can obtain the following equations: As lim t→∞ N h (t) = ∞, we only need to consider the following limit system: For system (12), there is only a equilibrium E 0 = (1, 0, 0, 0), and the equilibrium is globally asymptotically stable. Furthermore, let P (t) = (S h (t), E h (t), I h (t), C h (t)) be the positive equilibrium of system (10) at time t. We can obtain: Taking the equations of system (13) into N h (t) = S h (t) + E h (t) + I h (t) + C h (t), we can obtain: Hence, The indirect transmission rate of human population is saturating incidence or mass action incidence, and b h > d h . When system (2) is stable, the last four equations of system (1) will become the following equations: Before discussing the details, we make the following scaling transforms for system (14): We have, Due to lim t→∞ N h (t) = ∞, so we only need to consider the following limit system: Let P * = (s * h , e * h , i * h , c * h ) be the positive equilibrium of system (16), we have: For the stability of positive equilibrium P * of system (16), we have the following theorem.
Theorem 2.3. The positive equilibrium P * of system (16) is globally asymptotically stable.