1. Introduction

Brucellosis, also called Malta fever, is a highly contagious zoonosis caused by bacteria of the genus Brucella and is one of the most wide-spread infectious diseases in the world. Brucellosis is rarely fatal and the mortality of livestock due to the disease is practically zero. But brucellosis has important impacts on the livestock industry and the loss due to the fact that abortions can be very high if no control is applied. Human brucellosis is commonly caused by exposure to infected livestock or livestock products[32]. There is no recorded transmission of the infection among humans and they can very rarely infect animals. Human brucellosis remains the commonest zoonotic disease worldwide with more than half a million new cases annually.

The global epidemiology of the disease has drastically evolved over the past decade. Pappas et al[29] depict the global distribution of the disease before 2006 and draw a new global map of human brucellosis. The papers[4,32] describe the history and development of brucellosis in China till the early 2000s and briefly present the variation of epidemic situation, epidemiological characteristics, application of vaccines, and disease controls. After implementing comprehensive measures, great progress had been actually achieved in the prevention and control of brucellosis in China[32]. Brucellosis, however, still remains a serious public health issue. Since 1992, new foci of human brucellosis have emerged in China and the situation in certain regions is rapidly worsening, particularly in Inner Mongolia, one of the most important animal husbandry provinces of China. The brucellosis spreads all over Inner Mongolia and most of the confirmed cases are peasants and herdsmen[6]. The incidence of human brucellosis in Inner Mongolia increased from 3.42 per 100 000 in 2002 to 33.32 per 100 000 in 2006[15]. In recent years, Inner Mongolia has the largest number (40$\%$) of new annually confirmed cases of human brucellosis in China. Most of the cases in Inner Mongolia are attributed to Brucella melitensis (sheep brucella) (92$\%$) and then to Brucella abortus (5$\%$) due to the rapid development of sheep farming[30]. In Inner Mongolia, some of human brucellosis cases recover to susceptible under good treatment while the others become chronic cases[41]. In particular, in rural areas, due to the lower awareness of the disease and the limited medical resources, many infected cases miss the best time to treat and become chronic cases[30,32]. The prevention and control of brucellosis require substantial financial support. So, it is necessary to find a balance between the control effectiveness and the control resources and then to formulate an optimal control strategy with the very limited support from the government.

The use of mathematical formulations and models has emerged as a means to better describe, understand and predict the dynamical properties of brucellosis evolution and provide insight into the mechanisms that lead to those dynamics[1,2,3,5,10,16,14,26,44]. Most study is based on statistical models (static study) [2,3,10,26,34]. There are also several dynamical models proposed to explore the transmission of brucellosis[1,5,14,16,44].

Most existing models of brucellosis only consider transmission within a single livestock population and the transmission to humans is rarely considered[1,5,14]. According to the domestic situation of Mongolia, Zinsstag et al[44] present a livestock-to-human brucellosis transmission model, which considers brucellosis transmission within sheep and cattle populations and spreading to humans as additive components. Based on the model, the authors[44] estimate demographic (birth rate, mortality) and transmission (contact rates) parameters between livestock and livestock to humans for cost-effectiveness analysis of a nation-wide mass-vaccination programme. Recently, Hou et al[16] propose a brucellosis transmission model, which involves sheep population, human population and brucella in the environment. The results[16] show that vaccinating and disinfecting both young and adult sheep are appropriate and effective strategy to control brucellosis in Inner Mongolia of China.

As far as we know, available dynamic models of brucellosis transmission are developed without considering the limitedness of control resources. In this study, we develop and analyze a dynamic model, which considers brucellosis transmission within sheep population and the transmission to humans as additive components, and explore the optimal control of brucellosis in resources-limited setting. The paper is organized as follows. In Section 2, a new sheep-to-human brucellosis transmission model is formulated. The global dynamics of the model is explored in Section 3 with the help of asymptotic autonomous system theory and Lyapunov direct method. Section 4 formulates a multi-objective optimization problem in resources-limited setting and then, by the weighted sum method, transforms it into a scalar optimization problem with minimizing the cost for vaccination and health education and also the incidences of brucellosis both in sheep and in human. Section 5 characterizes the effects of sheep recruitment, vaccination of sheep, culling of infected sheep, and health education on the dynamics of brucellosis and on its optimal control. The paper ends with some conclusions and enlightening discussions.

---

2. Model formulation

Based on the transmission characteristics of brucellosis in Inner Mongolia of P. R. China, we consider both the sheep population and the human population. The sheep population is divided into three epidemiological classes, the susceptible ($S$), infectious ($I$), and vaccinated ($V$). The human population consists of a high risk subpopulation ($G_1$) and a low risk subpopulation ($G_2$) according to the infection risk. Each human subpopulation is partitioned into three classes, the susceptible, acute human brucellosis, and chronic human brucellosis, with population size denoted by $S_i$, $I_i$, $C_i$, $i=1,2$, respectively. The dynamic transfer of hosts is depicted in Fig. 1.

In sheep population, the susceptible admits effective vaccination at a constant rate $\theta$. The influx of susceptible comes from two sources, a constant recruitment $A$ and hosts that lost immunity $\delta V$. The natural elimination rate is assumed to be the same constant $m$ for all the hosts and the infectious hosts suffer an extra culling with constant rate $k$. The incidence term is of the bilinear mass-action form $\lambda SI$.

In human population, each subpopulation admits a constant recruitment $\Lambda_i$ and natural death at a constant rate $d_i$. $1/\gamma_{i}$ denotes the acute onset period of human brucellosis. $q_i$ is the fraction of acute human brucellosis turned into chronic cases. Both incidence terms of $G_1$ and $G_2$ are also of the bilinear type. The per capita transmission rate of $G_1$ is set to be $\beta$ and that of $G_2$ reads $\varepsilon\beta$ with $\varepsilon$ ($0<\varepsilon\ll1$) being the infection risk attenuation coefficient of $G_2$ compared to $G_1$.

The description, unit, default value and reference resource for each parameter are given in Table 1. The default value for each of the parameters are derived from the references[8,9,27,33,41,42] on brucellosis of P. R. China.

The transfer diagram in Fig. 1 leads to the following system of differential equations

$\left\{\begin{array}{ll} \begin{aligned}\label{SheepHumanModel} \frac{dS}{dt}&=A+\delta V-\lambda SI-mS-\theta S,\\ \frac{dI}{dt}&=\lambda SI-mI-k I,\\ \frac{dV}{dt}&=\theta S-\delta V-mV,\\ \frac{dS_1}{dt}&=\Lambda_1+(1-q_1)\gamma_1 I_1-\beta S_1I-d_1S_1,\\ \frac{dI_1}{dt}&=\beta S_1I-\gamma_1 I_1-d_1I_1,\\ \frac{dC_1}{dt}&=q_1\gamma_1 I_1-d_1C_1,\\ \frac{dS_2}{dt}&=\Lambda_2+(1-q_2)\gamma_2 I_2-\varepsilon\beta S_2I-d_2S_2,\\ \frac{dI_2}{dt}&=\varepsilon\beta S_2I-\gamma_2 I_2-d_2I_2,\\ \frac{dC_2}{dt}&=q_2\gamma_2 I_2-d_2C_2. \end{aligned} \end{array} \right.$

(1)

---

3. Global dynamics of (1)

In this section, the global dynamics of (1) will be explored. By the standard next generation method [37], the basic reproduction number of (1) reads

$R_{0}=\displaystyle\frac{\lambda(m+\delta)A}{m(m+k)(m+\delta+\theta)}.$

(2)

It is trivial to show that the following compact feasible region

$\label{feasible-region} \begin{split} \Gamma=\Big\{(S,I,V,S_1,I_1,C_1,S_2,I_2,C_2)\in\mathbb{R}_+^{9}\Big|& 0\leq S+I+V\leq \frac{A}{m},\\ 0\leq S_1+I_1+C_1\leq \frac{\Lambda_1}{d_1}, 0\leq S_2+I_2+C_2\leq \frac{\Lambda_2}{d_2} \Big\} \end{split}$

is positively invariant with respect to (1), and (1) always has the disease-free equilibrium

$E_{0}=(S^0,0,V^0,S_1^0,0,0,S_2^0,0,0),$

where

$S^0=\displaystyle\frac{(m+\delta)A}{m(m+\delta+\theta)},~~V^0=\displaystyle\frac{\theta A}{m(m+\delta+\theta)},~~S_i^0=\displaystyle\frac{\Lambda_i}{d_i},~~(i=1,2).$

(3)

Follow Theorem 2 of [37], it is not difficult to reach the following claim on the stability of $E_0$.

Theorem 3.1. The disease-free equilibrium $E_{0}$ is locally asymptotically stable if $R_0<1$, and is unstable if $R_0>1$.

Moreover, (1) has a unique endemic equilibrium

$E^*=(S^*,I^*,V^*,S_1^*,I_1^*,C_1^*,S_2^*,I_2^*,C_2^*)$

if and only if $R_0>1$, where

$\begin{split} S^*&=\frac{m+k}{\lambda}, ~~V^*=\frac{\theta}{m+\delta}S^*,\\ I^*&=\displaystyle\frac{m(m+\delta+\theta)}{\lambda(m+\delta)}\Big[\displaystyle\frac{A\lambda(m+\delta)}{m(m+\delta+\theta)(m+k)}-1\Big] =\displaystyle\frac{m(m+\delta+\theta)}{\lambda(m+\delta)}(R_0-1),\\ S_1^*&=\displaystyle\frac{\gamma_1I_1^*+d_1I_1^*}{\beta I^*},~~C_1^*=\displaystyle\frac{q_1\gamma_1I_1^*}{d_1},~~I_1^*=\displaystyle\frac{\beta\Lambda d_1I^*}{\beta q_1\gamma_1 I^*+\beta d_1 I^*+d_1(\gamma_1+d_1)}, \\S_2^*&=\displaystyle\frac{\gamma_2I_2^*+d_2I_2^*}{\varepsilon \beta I^*},~~C_2^*=\displaystyle\frac{q_2\gamma_2I_2^*}{d_2},~~I_2^*=\displaystyle\frac{\varepsilon\beta\Lambda d_2I^*}{\varepsilon\beta q_2\gamma_2 I^*+\varepsilon\beta d_2 I^*+d_2(\gamma_2+d_2)}. \end{split}$

(4)

Next, we explore the globally asymptotic stability of (1). The main approach involves the theory of asymptotic autonomous systems [23,36]. For the reader's convenience, we shall first summarize below a few concepts and results of asymptotic autonomous system from [23] that will be basic for the following discussion.

Consider the following differential equations

$\frac{dx}{dt}=f(t,x),$

(5)

$\frac{dy}{dt}=g(y).$

(6)

(5) is called asymptotically autonomous with limit equation (6) if

$f(t,x)\rightarrow g(x), t\rightarrow \infty,~~\text{locally uniformly in}~~x \in \mathbb{R}^n,$

i.e., for $x$ in any compact subset of $\mathbb{R}^n$. Assume that $f(t,x), g(x)$ are continuous functions and locally Lipschitz in $x$.

Lemma 3.2 ([23]). The $\omega$-limit set $\Omega$ of a forward bounded solution $x$ to (5) is non-empty, compact, and connected. Moreover, $\Omega$ attracts $x$, i.e.,

$dist(x(t),\Omega)\rightarrow 0, t\rightarrow \infty.$

Finally $\Omega$ is invariant under (6). In particular any point in $\Omega$ lies on a full orbit of (6) that is contained in $\Omega$.

Lemma 3.3 ([23]). Let $e$ be a locally asymptotically stable equilibrium of (6) and $\Omega$ the $\omega$-limit set of a forward bounded solution $x$ of (5). If $\Omega$ contains a point $y_0$ such that the solution $y$ of (6), with $y(0)=y_0$, converges to $e$ for $t\rightarrow \infty$, then $\Omega=\{e\}$, i.e., $x(t)\rightarrow e, t\rightarrow \infty$.

Now, consider the following sub-model of (1)

$\left\{\begin{array}{ll} \begin{aligned} \frac{dS}{dt}&=A+\delta V-\lambda SI-mS-\theta S,\\ \frac{dI}{dt}&=\lambda SI-mI-k I,\\ \frac{dV}{dt}&=\theta S-\delta V-mV. \end{aligned} \end{array} \right.$

(7)

Lemma 3.4. For system (7), we have the following claims:

(1) If $R_0\leq1$, then the disease-free equilibrium $(S^0,0,V^0)$ of (7) is globally asymptotically stable;

(2) If $R_0>1$, then the unique endemic equilibrium $(S^*,I^*,V^*)$ of (7) is globally asymptotically stable.

Proof. (1) Consider a Lyapunov function defined by

$L(S,I,V)=\Big(S-S^0-S^0\text{ln}\frac{S}{S^0}\Big)+\Big(V-V^0-V^0\text{ln}\frac{V}{V^0}\Big)+I.$

Differentiating $L(S,I,V)$ along the solutions of (7) and using the equilibrium equations $A=m S^0+\theta S^0-\delta V^0$ and $m V^0=\theta S^0-\delta V^0$, we obtain

$\begin{aligned} \frac{dL}{dt}=&{{(1-\frac{S^0}{S})\frac{dS}{dt}+(1-\frac{V^0}{V})\frac{dV}{dt}+\frac{dI}{dt}}}\\ =&(m+k)(R_0-1)I+A(2-\frac{S}{S^0}-\frac{S^0}{S})+\delta V^0(1+\frac{V}{V^0}-\frac{S}{S^0}-\frac{S^0V}{SV^0})\\&+\theta S^0(1+\frac{S}{S^0}-\frac{V}{V^0}-\frac{SV^0}{S^0V})\\ =&(m+k)(R_0-1)I+m S^0(2-\frac{S}{S^0}-\frac{S^0}{S})+\theta S^0(2-\frac{S}{S^0}-\frac{S^0}{S})\\ &-\delta V^0(2-\frac{S}{S^0}-\frac{S^0}{S})+\theta S^0(1+\frac{S}{S^0}-\frac{V}{V^0}-\frac{SV^0}{S^0V})\\&-\delta V^0(1+\frac{S}{S^0}-\frac{V}{V^0}-\frac{SV^0}{S^0V})+\delta V^0(1+\frac{S}{S^0}-\frac{V}{V^0}-\frac{SV^0}{S^0V})\\&+\delta V^0(1+\frac{V}{V^0}-\frac{S}{S^0}-\frac{S^0V}{SV^0})\\ =&(m+k)(R_0-1)I+m S^0(2-\frac{S}{S^0}-\frac{S^0}{S})+m V^0(3-\frac{S^0}{S}-\frac{V}{V^0}-\frac{SV^0}{S^0V})\\ &+\delta V^0(2-\frac{SV^0}{S^0V}-\frac{S^0V}{SV^0}). \end{aligned}$

If $R_0\leq1$, then $dL/dt\leq 0$ and $dL/dt=0$ if and only if $(S,I,V)=(S^0,0, V^0)$. Therefore, the largest compact invariant set in $\{(S, I, V)| dL/dt=0\}$ is the singleton $\{(S^0,0,V^0)\}$. The LaSalle's invariant principle [20] implies that $(S^0,0,V^0)$ is globally asymptotically stable when $R_0\leq1$.

(2) Consider the Lyapunov function

$L(S,I,V)=S-S^*-S^*\text{ln}\frac{S}{S^*}+I-I^*-I^*\text{ln}\frac{I}{I^*}+V-V^*-V^*\text{ln}\frac{V}{V^*}.$

Differentiating $L(S,I,V)$ along solutions of (7) produces

$\begin{aligned} \frac{dL}{dt}=&(1-\frac{S^*}{S})(A+\delta V-\lambda SI-mS-\theta S)+(1-\frac{I^*}{I})(\lambda SI-mI-k I)\\&+(1-\frac{V^*}{V})(\theta S-\delta V-mV)\\ =&A(2-\frac{S^*}{S}-\frac{S}{S^*})+\delta V^*(\frac{S}{S^*}-1)(\frac{S^* V}{S V^*}-1)+\theta S^*(\frac{V}{V^*}-1)(\frac{S V^*}{S^* V}-1)\\ =&2A+\delta V^*+\theta S^*-\frac{S}{S^*}(mS^*+\lambda S^*I^*)-mV-A\frac{S^*}{S}-\delta V\frac{S^*}{S}-\theta S\frac{V^*}{V}\\ =&(\lambda S^*I^*+mS^*)(2-\frac{S^*}{S}-\frac{S}{S^*})+\delta V^*(2-\frac{S^* V}{S V^*}-\frac{S V^*}{S^* V})\\&+3mV^*(3-\frac{S^*}{S}-\frac{V}{V^*}-\frac{S V^*}{S^* V})\leq 0. \end{aligned}$

Obviously, $dL/dt=0$ if and only if $(S,I,V)=(S^*,I^*,V^*)$. It implies that, when $R_0>1$, the largest compact invariant set of the sub-model (7) in $\{(S, I, V)|dL/dt=0\}$ is the singleton $\{(S^*,I^*,V^*)\}$. Therefore, the equilibrium $(S^*,I^*,V^*)$ of (7) is globally asymptotically stable.

Theorem 3.5. For system (1), we have the following conclusions:

(1) If $R_0\leq1$, then the disease-free equilibrium $E_{0}$ of (1) is globally asymptotically stable;

(2) If $R_0>1$, then the unique endemic equilibrium $E^*$ of (1) is globally asymptomatically stable.

Proof. $(1)$ Note that, the first three equations of (1) are independent of $S_i$, $C_i$ and $I_i$, $i=1,2$. Moreover, Lemma 3.4 asserts the globally asymptotic stability of the disease-free equilibrium for the first three equations of (1) when $R_0\leq1$.

Now, consider the rest equations of (1) except (7), that is,

$\left\{\begin{array}{lll} \begin{aligned} \frac{dS_1}{dt}&=\Lambda_1+(1-q_1)\gamma_1 I_1-\beta S_1I(t)-d_1S_1,\\ \frac{dI_1}{dt}&=\beta S_1I(t)-\gamma_1 I_1-d_1I_1,\\ \frac{dC_1}{dt}&=q_1\gamma_1 I_1-d_1C_1\\ \frac{dS_2}{dt}&=\Lambda_2+(1-q_2)\gamma_2 I_2-\varepsilon\beta S_2I(t)-d_2S_2,\\ \frac{dI_2}{dt}&=\varepsilon\beta S_2I(t)-\gamma_2 I_2-d_2I_2,\\ \frac{dC_2}{dt}&=q_2\gamma_2 I_2-d_2C_2. \end{aligned} \end{array} \right.$

(8)

When $R_0\leq1$, one has $\lim\limits_{{t\rightarrow\infty}}I(t)=0$. Then, the limiting system of (8) is

$\left\{\begin{array}{ll} \begin{aligned} {{\frac{dS_1}{dt}}}&{{=\Lambda_1+(1-q_1)\gamma_1 I_1-d_1S_1,}}\\ {{\frac{dI_1}{dt}}}&{{=-\gamma_1 I_1-d_1I_1,}}\\ {{\frac{dC_1}{dt}}}&{{=q_1\gamma_1 I_1-d_1C_1,}}\\ {{\frac{dS_2}{dt}}}&{{=\Lambda_2+(1-q_2)\gamma_2 I_2-d_2S_2,}}\\ {{\frac{dI_2}{dt}}}&{{=-\gamma_2 I_2-d_2I_2,}}\\ {{\frac{dC_2}{dt}}}&{{=q_2\gamma_2 I_2-d_2C_2}}. \end{aligned} \end{array} \right.$

(9)

Obviously, the linear system (9) admits a unique equilibrium $(\Lambda_1/d_1,0,0,\Lambda_2/d_2,0,0)$ being globally asymptotically stable.}} From the theory of asymptotic autonomous systems, it follows that the disease-free equilibrium of (8) $(\Lambda_1/d_1,0,0,\Lambda_2/d_2,0,0)$ is globally asymptotically stable.

(2) When $R_0>1$, the limiting system of (8) is

$\left\{\begin{array}{lll} \begin{aligned} \frac{dS_1}{dt}&=\Lambda_1+(1-q_1)\gamma_1 I_1-\beta S_1I^*-d_1S_1,\\ \frac{dI_1}{dt}&=\beta S_1I^*-\gamma_1 I_1-d_1I_1,\\ \frac{dC_1}{dt}&=q_1\gamma_1 I_1-d_1C_1\\ \frac{dS_2}{dt}&=\Lambda_2+(1-q_2)\gamma_2 I_2-\varepsilon\beta S_2I^*-d_2S_2,\\ \frac{dI_2}{dt}&=\varepsilon\beta S_2I^*-\gamma_2 I_2-d_2I_2,\\ \frac{dC_2}{dt}&=q_2\gamma_2 I_2-d_2C_2. \end{aligned} \end{array} \right.$

(10)

It is not difficult to show that the linear system (10) has an equilibrium $(S_1^*,I_1^*,C_1^*,S_2^*,I_2^*,C_2^*)$ being globally asymptotically stable. By the theory of asymptotic autonomous systems, the equilibrium $(S_1^*,I_1^*,C_1^*,S_2^*,I_2^*,C_2^*)$ of (8) is globally asymptotically stable.

---

4. Optimal control strategy

The brucellosis asserts heavy burden to human health and local economy. The government consequently has to adopt a series of measures to control the brucellosis. In theory, by the global threshold dynamics of (1), from the explicit expressions of $R_0$ in (2) and $I^*, I_1^*, I_2^*$ in (4), it follows that $R_0, I^*, I_1^*, I_2^*$ decrease with the decreasing of $\lambda, \delta, A$ and with the increasing of $k, \theta, m$. That is to say, in order to control the brucellosis, it is reasonable to reduce the transmission between sheep (e.g., improve the living conditions of sheep), effectively inoculate more sheep with vaccine having long immunization period, reduce the size of sheep population, remove (cull or naturally eliminate) more infectious sheep. In the long run, all those measures benefit the control of brucellosis and should be policy suggestions. It may be so in theory, but how will it work in practice?

In reality, the financial support from the government is very limited and can not meet the real requirement of the brucellosis control. The testing of sheep brucella and the culling of infected sheep can not be sustainedly and effectively implemented in large scale[41]. Some moderate measures have been carried out to control the brucellosis such as vaccination program in sheep stock and health education in human population aiming to raise the consciousness among publics to protect against brucellosis by influencing self examination behavior.

Taking into account the serious situation of the brucellosis in Inner Mongolia and the limited financial support from the government, one has to find a balance between the loss induced by brucellosis and the cost of control effort. It is very important and reasonable to achieve the best control effect with given limited financial budget.

Optimization has been playing a key role in the design, planning and operation of infectious disease control and related processes. It is reasonable and challenging to develop optimal strategy for more effective control or treatment options of infectious diseases in resource-limited settings. Brucellosis infections and their control are a world-wide challenge due to limited available resources.

In the following, we explore the optimal control strategy for brucellosis by considering both the control effectiveness and resource limitation. According to the reality of brucellosis control in Inner Mongolia, we still take vaccination in sheep stock and health education in human population as control measures in our optimal control model. Then the vaccination rate $\theta$, in (1), is set to be time-dependent $\theta(t)$ describing the time-varying vaccination strategy. The functions $\varphi_i(t)$ measure the effort of health education in subpopulation $G_i (i=1,2)$, respectively. It means that $\varphi_i(t)S_i(t) (i=1,2)$ are of susceptible who have acquired effective brucellosis prevention measures (i.e., only $\beta(1-\varphi_1(t))S_1(t)I(t)$ and $\varepsilon\beta(1-\varphi_2(t))S_2(t)I(t)$ susceptible are infected brucellosis in high and low risk subpopulation per unit time, respectively). Since it seems impossible to vaccinate all the susceptible sheep and to educate the entire human population due to the limited resources, it is realistic to assume that there exist $\bar{\theta}, \bar{\varphi_i}\in (0,1)$ such that $0\leq \theta(t)\leq \bar{\theta}<1$ and $0\leq \varphi_i(t)\leq \bar{\varphi_i}<1, i=1,2$. In the fixed time interval $[0,T]$ with $T>0$, the feasible decision space or constraint set reads

$\begin{split} \mathbf{U}=\Big\{(\theta(t), \varphi_1(t), \varphi_2(t))\Big|&0\leq \theta(t)\leq \bar{\theta},~0\leq \varphi_i(t)\leq \bar{\varphi_i},~0\leq t \leq T,~ \\&\theta(t) \text{and}~\varphi_i(t) \text{are Lebesgue measurable}, i=1,2\Big\}. \end{split}$

(11)

In most real-life optimization scenarios and designs of the control of brucellosis, multiple objectives under consideration arise naturally and often conflict with each other. Anyone actually prefers to achieve best control by using less effort or cost, i.e., minimize cost, maximize performance, maximize reliability, etc. It seems very difficult but realistic.

Here our key objectives are to minimize the cost for vaccination and health education and to minimize the economic loss caused by culling of infected sheep and treatment of infected human cases. These costs during $[0, T]$ read

$\begin{array}{ll} \displaystyle\int_0^TB_0S(t)\theta(t)dt,&\displaystyle\int_0^T\sum\limits_{i=1}^2B_iS_i(t)\varphi_i(t)dt,\\ \displaystyle\int_0^T D_0kI(t)dt,&\displaystyle\int_0^T\Big(D_1(I_1(t)+I_2(t))+D_2(C_1(t)+C_2(t))\Big)dt, \end{array}$

(12)

respectively, where $B_0$ is the cost of brucellosis vaccine for each sheep, $B_i$ is the cost of health education per unit time for subpopulation $G_i(i=1,2)$, $D_0$ denotes the economic loss of culling each infected sheep, $D_1$ and $D_2$ denote the cost of medical resources for each acute and chronic human brucellosis case respectively.

Now, we are facing with a multi-objective optimization problem

$\begin{array}{ll} \min\displaystyle\int_0^T D_0kI(t)dt,&\min\displaystyle\int_0^T\Big(D_1(I_1(t)+I_2(t))+D_2(C_1(t)+C_2(t))\Big)dt,\\ \min\displaystyle\int_0^TB_0S(t)\theta(t)dt,&\min\displaystyle\int_0^T\sum\limits_{i=1}^2B_iS_i(t)\varphi_i(t)dt, \end{array}$

(13)

subject to

$\begin{split} &\displaystyle\frac{dS}{dt}=A+\delta V-\lambda SI-mS-\theta(t) S,\\ &\displaystyle\frac{dI}{dt}=\lambda SI-mI-k I,\\ &\displaystyle\frac{dV}{dt}=\theta(t) S-\delta V-mV,\\ &\displaystyle\frac{dS_1}{dt}=\Lambda_1+(1-q_1)\gamma_1 I_1-\beta(1-\varphi_1(t)) S_1I-d_1S_1,\\ &\displaystyle\frac{dI_1}{dt}=\beta(1-\varphi_1(t)) S_1I-\gamma_1 I_1-d_1I_1,\\ &\displaystyle\frac{dC_1}{dt}=q_1\gamma_1 I_1-d_1C_1,\\ &\displaystyle\frac{dS_2}{dt}=\Lambda_2+(1-q_2)\gamma_2 I_2-\varepsilon\beta(1-\varphi_2(t)) S_2I-d_2S_2,\\ &\displaystyle\frac{dI_2}{dt}=\varepsilon\beta(1-\varphi_2(t)) S_2I-\gamma_2 I_2-d_2I_2,\\ &\displaystyle\frac{dC_2}{dt}=q_2\gamma_2 I_2-d_2C_2,\\ &(\theta(t), \varphi_1(t), \varphi_2(t))\in \mathbf{U}, \end{split}$

(14)

with initial conditions

$\begin{split} &0\leq S(0)\leq M, 0\leq I(0)\leq M, 0\leq V(0)\leq M, \\&0\leq S_i(0)\leq M, 0\leq I_i(0)\leq M, 0\leq C_i(0)\leq M,~~i=1,2, \end{split}$

(15)

where $M$ is a positive constant.

In term of methodology, there are several different approaches to solve general continuous multi-objective optimization problems[24]. One of the most common approaches is in general known as the weighted sum or scalarization method, which minimizes a positively weighted convex sum of the objectives by combining its multiple objectives into one single-objective scalar composite function with equal or weighted treatment.

By the weighted sum method, we transformed the vector-valued optimization problem (13)-(15) into a scalar quadratic optimization problem with a unique objective function of the following form

$\begin{array}{ll} J\Big(\theta(t), \varphi_1(t), \varphi_2(t)\Big)=&\displaystyle\int_0^T\Bigg[\omega_1\Big(D_0kI(t)\Big)^2+\omega_2\Big(\Big(D_1(I_1(t)+I_2(t))\Big)^2\\ &+\Big(D_2(C_1(t)+C_2(t))\Big)^2\Big)+\omega_3\Big(B_0S(t)\theta(t)\Big)^2\\&+\omega_4\sum_{i=1}^2\Big(B_iS_i(t)\varphi_i(t)\Big)^2\Bigg]dt, \end{array}$

(16)

subjected to (14), (15), and $\sum\limits_{i=1}^4\omega_i=1$, $\omega_i\in [0,1]$. Here, the objective function $J\Big(\theta(t), \varphi_1(t), \varphi_2(t)\Big)$ measures the total economic loss caused by the brucellosis.

To facilitate the discussions below, let $X=(S,I,V,S_1,I_1,C_1,S_2,I_2,C_2)'$, $u=(\theta(t),\varphi_1(t),\varphi_2(t))'$, $F(t,X,u)=(f_1,f_2,...,f_9)'$, where the prime represents the transpose of a row vector and $f_i$ is the $i$-th right-hand side function of (14). Then the above single objective optimization problem rewrites

$\begin{split} \min\limits_{(\theta(t), \varphi_1(t), \varphi_2(t))\in \mathbf{U}} J(\theta(t),\varphi_1(t),\varphi_2(t)), \end{split}$

(17)

subject to

$\displaystyle\frac{dX}{dt}=F(t,X,u),~~X(0)\geq 0,$

(18)

where $\|X(0)\|\leq 3M$ and $\omega_i\in [0,1]$ such that $\sum\limits_{i=1}^4\omega_i=1$.

In the optimal control theory, after formulating the model and the corresponding objective functionals appropriate to the scenarios, one has to explore several basic problems such as proving the existence of an optimal control and characterizing the optimal control. In order to prove the existence of an optimal control, we introduce the following conclusion ([12], Theorem 4.1 and Corollary 4.1).

Lemma 4.1. Consider the optimal control problem

$\min\limits_{u\in \mathbf{D}}J(u)=\min\limits_{u\in \mathbf{D}}\int_{t_0}^{t_1}L(t,x,u)dt+\phi(e)$

subject to $\displaystyle\frac{dx}{dt}=f(t,x,u)$ with $e=(t_0,t_1,x(t_0),x(t_1))\in S\subset\mathbb{R}^{2n+2}$, where $L(t,x,u):\mathbb{R}^1 \times \mathbb{R}^n \times \mathbb{R}^m\rightarrow\mathbb{R}^1$, $\phi(e):\mathbb{R}^1\times \mathbb{R}^1\times \mathbb{R}^n\times \mathbb{R}^n\rightarrow\mathbb{R}^1$, and $f(t,x,u):\mathbb{R}^1\times \mathbb{R}^n\times \mathbb{R}^m\rightarrow\mathbb{R}^n$ are continuous and $f(t,x,u)$ satisfies Lipschitz property with respect to $x$, $\mathbf{D}=\{u(t)=(u_1(t),u_2(t),...,u_m(t)): u_i(t) is~Lebesgue~ measurable~on~[t_0,t_1], i=1,2,...,m\}.$ Assume that, $(i)$ $f(t,x,u)$ is $C^{1}$, the partial derivatives of $f(t,x,u)$ and $f(t,0,0)$ are bounded, $(ii)$ The set of control $u$ and corresponding state variable $x$ is non-empty, $(iii)$ $\mathbf{D}$ is closed, $(iv)$ $S$ is compact and $\phi$ is continuous on $S$, $(v)$ There exist constants $C_1>0$, $C_2$, and $\gamma>1$ such that $L(t,x,u)\geq C_1|u|^{\gamma}-C_2$, (vi) $\mathbf{D}^+=\Big\{(n,m)\Big|\exists v\in\mathbf{D}, m=f(t,x,v),~n\geq L(t,x,v)\Big\}$ is convex for $(t,x)$. Then there exists an optimal control $u^*\in\mathbf{D}$.

Theorem 4.2. There exists an optimal pair $(u^*,X^*)$ such that

$J(\theta^*(t),\varphi_1^*(t),\varphi_2^*(t))=\min\limits_{(\theta, \varphi_1, \varphi_2)\in \mathbf{U}} J(\theta(t),\varphi_1(t),\varphi_2(t))$

subject to (18), where $u^*=(\theta^*(t), \varphi_1^*(t),\varphi_2^*(t))$ is the optimal control, $X^*=(S^*(t)$, $I^*(t)$, $V^*(t)$, $S_1^*(t)$, $I_1^*(t)$, $C_1^*(t)$, $S_2^*(t)$, $I_2^*(t)$, $C_2^*(t))$ is the state solution corresponding to the optimal control.

Proof. It suffices to show that the conditions of Lemma 4.1 are satisfied with $\phi(0,T,X(0),X(T))=0$ since (17)-(18) is a special case of the optimal control problem in Lemma 4.1 with $\phi(e)\equiv 0$.

The claim $(i)$ of Lemma 4.1 follows from the properties of $F(t,X,u)$ and the state variables of (18). The existence of solutions of (18) with $u\in \mathbf{U}$ is guaranteed by Theorem 9.2.1 in [19], which implies $(ii)$ of Lemma 4.1. The constraint set $\mathbf{U}$ is closed and convex since $u\in\mathbf{U}$ is bounded and Lebesgue measurable (i.e., $L^\infty$ bound hold). $(iv)$ of Lemma 4.1 holds due to $\phi(0,T,X(0),X(T))=0$ and the set of $X(0)$ is bounded and closed.

Let $\Phi(t,X,v)=\omega_1\Big(D_0kI(t)\Big)^2+\omega_2\Big(\Big(D_1(I_1(t)+I_2(t))\Big)^2$$+\Big(D_2(C_1(t)+C_2(t))\Big)^2\Big) +\omega_3\Big(B_0S(t)\theta(t)\Big)^2+\omega_4\sum_{i=1}^2\Big(B_iS_i(t)\varphi_i(t)\Big)^2$. It is not difficult to verify that there exist $C_1>0$, $C_2$, and $\gamma>1$ such that

$\Phi(t,X,v)\geq C_1(|\theta|^2+|\varphi_1|^2+|\varphi_2|^2)^{\frac{\gamma}{2}}-C_2$

since $\theta(t)$, $\varphi_1(t)$, and $\varphi_2(t)$ are bounded.

In order to verify $(vi)$ of Lemma 4.1, set

$\begin{split} \mathbf{U}^+=\Big\{(n,m)\Big|&\exists v\in\mathbf{U}, m=F(t,X,v),~n\geq \Phi(t,X,v)\Big\}. \end{split}$

For any $(n_i,m_i)\in\mathbf{U}^+$, $i=1,2$, i.e., there exist $v_1\in\mathbf{U}$ and $v_2\in\mathbf{U}$ such that

$\begin{split} m_1=F(t,X,v_1),~~n_1\geq\Phi(t,X,v_1),\\ m_2=F(t,X,v_2),~~n_2\geq\Phi(t,X,v_2). \end{split}$

Since $F(t,X,u)$ is linear in $u$, for any $s\in[0,1]$,

$\begin{split} sm_1+(1-s)m_2=sF(t,X,v_1)+(1-s)F(t,X,v_2)=F(t,X,sv_1+(1-s)v_2). \end{split}$

Note that the control set $\mathbf{U}$ is convex. Then, there exists a $u\in \mathbf{U}$ such that $u=sv_1+(1-s)v_2$ and

$sm_1+(1-s)m_2=F(t,X,sv_1+(1-s)v_2)=F(t,X,u)$

for $s\in[0,1]$. On the other hand, since $\Phi(t,X,u)$ is convex in $u$, one has

$\begin{split} sn_1+(1-s)n_2&\geq s\Phi(t,X,v_1)+(1-s)\Phi(t,X,v_2)\\&\geq \Phi(t,X,sv_1+(1-s)v_2)=\Phi(t,X,u). \end{split}$

Then, for $s\in[0,1]$, $s(n_1,m_1)+(1-s)(n_2,m_2)=\Big(sn_1+(1-s)n_2,sm_1+(1-s)m_2\Big)\in\mathbf{U}^+.$ That is, $\mathbf{U}^+$ is convex. The proof is complete.

Next, we characterize the optimal control by deriving necessary conditions for the optimal control with the help of Pontryagin's Maximum Principle.

Theorem 4.3. The optimal control of (17) is characterized by

$\begin{split} \theta^*(t)&=min\Big\{\bar{\theta},max\Big\{\displaystyle\frac{\lambda_1-\lambda_3}{2\omega_3B_0^2S^*(t)},0\Big\}\Big\},\\ \varphi_1^*(t)&=min\Big\{\bar{\varphi_1},max\Big\{\displaystyle\frac{(\lambda_5-\lambda_4)\beta I^*(t)}{2\omega_4B_1^2S^*_1(t)},0\Big\}\Big\},\\ \varphi_2^*(t)&=min\Big\{\bar{\varphi_2},max\Big\{\displaystyle\frac{(\lambda_8-\lambda_7)\varepsilon \beta I^*(t)}{2\omega_4 B_2^2S_2^*(t)},0\Big\}\Big\},\\ \end{split}$

(19)

where $\lambda_i, i=1,2,...,9$ are the adjoint variables satisfying the adjoint equations (20) and the transversality conditions (21).

Proof. In order to facilitate the discussion below, we kill the stars in the optimal state $X^*(t)=(S^*(t), I^*(t), V^*(t), S^*_1(t), I^*_1(t), C^*_1(t), S^*_2(t), I^*_2(t), C^*_2(t))'$ corresponding to the optimal control $u^*=(\theta^*(t),\varphi_1^*(t),\varphi_2^*(t))$ and simply denote it by $(S(t), I(t), V(t), S_1(t), I_1(t), C_1(t), S_2(t), I_2(t), C_2(t))'$.

By Pontryagin's Maximum Principle, as in [7,28,31], to find the optimal control of (17)-(18) is equivalent to minimize the following Hamiltonian

$\begin{split} H=&\omega_1\Big(D_0kI(t)\Big)^2+\omega_2\Big(\Big(D_1(I_1(t)+I_2(t))\Big)^2+\Big(D_2(C_1(t)+C_2(t))\Big)^2\Big)\\ &+\omega_3\Big(B_0S(t)\theta(t)\Big)^2+\omega_4\sum\limits_{i=1}^2\Big(B_iS_i(t)\varphi_i(t)\Big)^2+\sum\limits_{i=1}^{9}\lambda_if_i, \end{split}$

where $\lambda_i, i=1,2,...,9,$ are the adjoint variables satisfying the following adjoint equations

$\begin{split} \lambda_1'=-\frac{\partial H}{\partial S}=&-2\omega_3B_0^2\theta^2(t)S(t)+\lambda_1 \lambda I+\lambda_1 m+\lambda_1 \theta(t)-\lambda_2 \lambda I-\lambda_3 \theta(t),\\ \lambda_2'=-\frac{\partial H}{\partial I}=&-2\omega_1D_0^2k^2I(t)+\lambda_1 \lambda S-\lambda_2(\lambda S-m-k)+\lambda_4 \beta(1-\varphi_1(t))S_1\\ &-\lambda_5 \beta(1-\varphi_1(t))S_1+\lambda_7 \varepsilon\beta(1-\varphi_2(t))S_2-\lambda_8 \varepsilon\beta(1-\varphi_2(t))S_2,\\ \lambda_3'=-\frac{\partial H}{\partial V}=&-\lambda_1 \delta +\lambda_3(\delta+m ),\\ \lambda_4'=-\frac{\partial H}{\partial S_1}=&-2\omega_4B_1^2\varphi_1^2(t)S_1(t)+\lambda_4(\beta(1-\varphi_1(t))I+d_1)-\lambda_5\beta(1-\varphi_1(t))I,\\ \lambda_5'=-\frac{\partial H}{\partial I_1}=&-2\omega_2D_1^2(I_1(t)+I_2(t))-\lambda_4(1-q_1)\gamma_1+\lambda_5(\gamma_1+d_1 )-\lambda_6q_1\gamma_1,\\ \lambda_6'=-\frac{\partial H}{\partial C_1}=&-2\omega_2D_2^2(C_1(t)+C_2(t))+\lambda_6 d_1,\\ \lambda_7'=-\frac{\partial H}{\partial S_2}=&-2\omega_4B_2^2\varphi_2^2(t)S_2(t)+\lambda_7 (\varepsilon\beta(1-\varphi_2(t))I+d_2)-\lambda_8\varepsilon\beta(1-\varphi_2(t))I,\\ \lambda_8'=-\frac{\partial H}{\partial I_2}=&-2\omega_2D_1^2(I_1(t)+I_2(t))-\lambda_7 (1-q_2)\gamma_2+\lambda_8(\gamma_2+d_2)-\lambda_9q_2\gamma_2,\\ \lambda_9'=-\frac{\partial H}{\partial C_2}=&-2\omega_2D_2^2(C_1(t)+C_2(t))+\lambda_9 d_2 \end{split}$

(20)

and the transversality conditions

$\begin{split} \lambda_i(T)=0,~~i=1,2,...,9. \end{split}$

(21)

The optimality conditions

$\displaystyle\frac{\partial H}{\partial \theta(t)}\Bigg|_{(\theta^*,\varphi^*_1(t),\varphi^*_2(t))}=0,~~\displaystyle\frac{\partial H}{\partial \varphi_i(t)}\Bigg|_{(\theta^*,\varphi^*_1(t),\varphi^*_2(t))}=0$

lead to

$\theta^*(t)=\displaystyle\frac{\lambda_1-\lambda_3}{2\omega_3B_0^2S^*(t)},~~\varphi_1^*(t)=\displaystyle\frac{(\lambda_5-\lambda_4)\beta I^*(t)}{2\omega_4B_1^2S^*_1(t)},~~\varphi_2^*(t)=\displaystyle\frac{(\lambda_8-\lambda_7)\varepsilon \beta I^*(t)}{2\omega_4 B_2^2S_2^*(t)}.\\ $

Moreover, taking into account the fact that $(\theta^*,\varphi^*_1(t),\varphi^*_2(t))\in\mathbf{U}$, i.e., using the lower and upper bounds of $u^*$, we find that the optimal control $u^*=(\theta^*(t),\varphi_1^*(t),\varphi_2^*(t))$ is well characterized by (19).

---

5. Sensitivity analyses and numerical simulations

In this section, we first use (1) to simulate the reported sheep brucellosis and human brucellosis data of Inner Mongolia from 2001 to 2010 and to predict the trends of the brucellosis. Then, we perform sensitivity analysis of some key parameters, numerically simulate the optimal control, and seek for some effective control and prevention measures.

The data of sheep population is extracted from China Animal Husbandry Yearbook [8] (see Table 2). The data, concerning human brucellosis from 2001 to 2014, are adopted from China Health Statistics Yearbook[25] (see Table 3). The data of sheep brucellosis can not be acquired easily since there are few published studies on the population dynamics of sheep, and we calculate the number of brucellosis sheep in Inner Mongolia (see Table 2) from data of sheep population in Table 2 and the annual seroprevalence of sheep brucella in key monitoring regions of Inner Mongolia[11,13,22,30,38,39,43].

In Inner Mongolia, the total population was about 24 million in the past decade with the birth rate less than $1\%$, and the rural population accounted for nearly half of the total population[25]. Since the vast majority of human brucellosis cases in Inner Mongolia are herdsmen, farmers, veterinarians, slaughters, and milkers, who live in rural region. So, we view the rural population in Inner Mongolia as the high risk subpopulation $G_1$ while the urban population in Inner Mongolia is considered as the low risk subpopulation $G_2$. Recent studies (e.g., [30]) revealed that, before 1980s, $94\%$ of human brucellosis cases in Inner Mongolia were from high risk subpopulation $G_1$, and, after 1990s, this percentage decreased to $85\%$. Therefore, for the current brucellosis outbreak in Inner Mongolia, it is reasonable to assume that $85\%$ of new confirmed human brucellosis are from the high risk subpopulation $G_1$ and the rest $15\%$ are from the low risk subpopulation $G_2$.

In order to carry out numerical simulations, most parameters of (1) are obtained from literature or assumed on the basis of common sense(see Table 1). The transmission rates of brucellosis (i.e., $\lambda$, $\beta$ and $\varepsilon$, see Table 1) are estimated by the nonlinear least-square method that produced the best fit between the model outputs and the data (see Fig. 2).

With the help of (1), we fit it with the data from 2001 to 2011 and predict the dynamics of both sheep brucellosis infection and human brucellosis infection in Inner Mongolia (Fig. 2). The numerical simulations show that the model (1) with reasonable parameter values provides a good match to the reported data. From 2001 to 2011, the brucellosis sheep and human brucellosis cases kept increasing. From 2012 to 2015, due to the government's strong intervention on the transmission of brucellosis (e.g., a special fund was set to cull the infected sheep) [17,35,40], the size of brucellosis sheep and human brucellosis cases had a significant decline. However, the epidemic situation will become more serious in the coming decades and reach a peak at about 2040 if the government of Inner Mongolia can not sustainedly provide enough support to cull the infected sheep (see Fig. 2).

In order to better understand the mechanism and the control of brucellosis, we characterize the influence of some key parameters or factors in (1) by sensitivity analysis. Fig. 3-6 provide sensitivity analysis for (1) in terms of $\theta$, $A$, $k$, and $\beta$, respectively.

From Fig. 3, we find that the vaccination of sheep can effectively control the brucellosis in some extent but can not eradicate it from either the sheep or human. Fig. 4 shows the sensitivity of (1) to $A$. It reveals that the population size of sheep has potential effect on the control of brucellosis. In fact, the number of sheep in Inner Mongolia has increased about $2.042\times10^7$ from 2001 to 2007. Reducing the breading scale of sheep and keeping the population size of sheep in a reasonable range or dividing the sheep into small sub-populations will facilitate the control of brucellosis.

Culling the infected sheep is another effective measure to control the brucellosis. Fig. 5 depicts that increasing the culling rate $k$ of infected sheep could effectively reduce the incidence of brucellosis in both sheep and human and can even eradicate the brucellosis if the culling rate exceeds some certain threshold level. Fig. 6 shows that reducing the transmission rate $\beta$ between the sheep and the high risk subpopulation $G_1$ facilitates the control of brucellosis but can not eliminate it. One typical way commonly used to effectively reduce $\beta$ is to carry out health and prevention education, which can enhance people' awareness and knowledge about brucellosis and the emergency measure and treatment to decrease the rate of clinical outbreak. From Fig. 6, we can see that the epidemic of brucellosis can be relieved, but cannot be eradicated with the current prevention and control measures.

According to the above sensitivity analysis and numerical experiments, by comparing those figures, the recruitment of susceptible sheep and the culling of infected sheep are more sensitive than vaccination and health education in the current brucellosis outbreak in Inner Mongolia. In conclusion, controlling the population of sheep, reducing the birth rate of sheep, increasing the vaccination rate of sheep, improving the culling management of infected sheep, enhancing the people's awareness of brucellosis, and any combination of these measures are effective to control brucellosis. In practice, some policy involving comprehensive measures must be brought into operation.

By Forward-Backward sweep method [21], we further expound the optimal control by numerical experiments. According to the current economic status in Inner Mongolia, the cost of culling each infected sheep is about 1000 CNY ($D_0=1$), the treatment cost of each acute and chronic human case is about 5000 CNY ($D_1=5$) and 3000 CNY ($D_2=3$) per year respectively, the vaccination of each susceptible sheep needs about 200 CNY ($B_0=0.2$), and the cost of health education for each person of high and low risk subpopulation are 300 CNY ($B_1=0.3$) and 200 CNY ($B_2=0.2$), respectively. Fig. 7 simulates the optimal control functions $\theta^*(t)$, $\varphi_1^*(t)$, $\varphi_2^*(t)$, and the corresponding states variables $I(t)^*$, $I_1^*(t)$, $I_2^*(t)$. The numeral simulations indicate that, with the current prevention and control strategy such as keeping sheep population in large size ($A=3300$) and culling rate of infected sheep at low level ($k=0.1$), putting majority of control resources into vaccination program and education plan, etc., the epidemic of brucellosis can be relieved in the short run, but can not be better controlled and eradicated in the long run. Even though $85\%$ of sheep population are effectively vaccinated and $90\%$ of human population admit a good prevention and health education (i.e., maximal supply of control resources), the number of brucellosis sheep still has a significant increase, the human brucellosis infection cases would continue decreasing in the next couple years but can not be eliminated(see Fig. 7).

The optimal control strategy explored in Section 4 sensitively depends on the selection of weighted coefficients. Since the principal aim is to reduce the number of brucellosis human and sheep cases, it is reasonable to assume that both $\omega_1$ and $\omega_2$ are greater than $\omega_3$ and $\omega_4$. For simplicity, we also could assume that the vaccination programme and the health education plan are equally important, i.e., $\omega_3=\omega_4$ when searching for an optimal control strategy. Fig. 8 illustrates the optimal control strategy with varying weighted coefficients. From the numerical simulation, we observe that the weighted coefficients have significant effect on the optimal control strategy. When $\omega_1$ increases from $0$ to $0.5$, the key minimizing targets transfer from vaccination and health education to culling of brucellosis sheep and treatment of human cases, and then the optimal control gradually increases from $\theta^*(t)=\varphi_1^*(t)=\varphi_2^*(t)=0$ to $\theta^*(t)=\overline{\theta}$, $\varphi_1^*(t)=\overline{\varphi}_1$ and $\varphi_2^*(t)=\overline{\varphi}_2$. Fig. 9 illustrates the varying of the total economic loss (i.e., the sum of the four integrals in (12)) and prevalence rates with respect to the weighted coefficients. In real application, the policy-maker should confirm some reasonable weighted coefficients according to the financial budget for brucellosis (e.g., TELoss$^*$ in Fig. 9) and the prevalence rates in both sheep and human population (e.g., $P_s^*,P_{h1}^*,P_{h2}^*$ in Fig. 9)).

Simulations in Fig. 10-11 reveal that the culling of infected sheep and the recruitment of sheep play important roles on the optimal control of brucellosis. Fig. 10 shows that, when the culling rate $k$ increases from $k=0.15$ to $k=0.7$, the brucellosis in both sheep and human is effectively controlled. Fig. 11 illustrates the role of the recruitment of sheep. When the constant influx rate of sheep decreases from $A=3300$ to $A=1600$, the brucellosis in both sheep and human can be better controlled.

---

6. Discussion and conclusions

The transmission of brucellosis has been a growing public concern in China, particularly in Inner Mongolia. In this study, we formulate a multigroup epidemiological model to study the transmission dynamics and explore control strategies for the brucellosis in Inner Mongolia. Our theoretical and numerical findings indicate that the brucellosis gradually increases in the next decades and reaches a peak at about 2030. Vaccination and health-prevention education, which are the main currently working control measures, positively affect and facilitate the disease control (see Fig. 3 and Fig. 6), but can not eliminate the brucellosis in Inner Mongolia (see Fig. 7). Sensitivity analyses in Fig. 4 and Fig. 5 show that the population size of sheep and the culling rate of brucellosis sheep are two important factors triggering the serious epidemic situation. The brucellosis in both sheep and human can be well controlled when the breeding size of sheep is reduced or the culling rate of infectious sheep is enhanced. Fig. 10 and Fig. 11 provide a series of effective and time-variable control strategies when the culling rate of brucellosis sheep is increased and when the population size of sheep is decreased, respectively. The numerical solutions of the optimal controls illustrated in Fig. 10 and Fig. 11 provides practical guidance for policymakers to make effective control policy and to implement it in resources-limited setting.

The study reveals that the brucellosis in Inner Mongolia of China is hardly eliminated if the policymakers still pay most attentions to the vaccination in sheep stock and the health-prevention education in human population. Our studies suggest that the government and the policymakers must take a new look at the current working control strategies. A reasonable and effective control strategy must involve comprehensive control measures based on the optimal control study.

Finally, we would like to point out that, even if a costly policy can prevent the disease spread, the prevention policy might not be well executed in resource-limited setting. A fundamental challenge in many control problems in epidemics is to find the right balance amongst several objectives. Although the weighted sum method is a typical method or technique to solve the multi-objective optimization, the problem lies in the correct selection of the weights to characterize the decision-makers preferences. In practice, it can be very difficult to precisely and accurately select these weights, even for someone very familiar with the problem domain. In addition, perturbations in the weights can lead to different solutions (see (19)). For this reason and others, decision-makers often prefer a set of promising solutions given the multiple objectives. In addition, recently, some evolutionary approaches have been proposed to deal with multi-objective optimization problems. Genetic algorithms, being a population based approach, is a well suited heuristic approach, see[18] for an overview and tutorial of multiple-objective optimization methods using genetic algorithms (GA). It is interesting but challenging to investigate the optimal control problem in this study by GA and compare it with the weighted sum method. It is left for our future work.

---

Acknowledgments

This work is partially supported by the National Natural Science Foundation of P. R. China (No. 11671072,11271065,11426045), the Research Fund for the Doctoral Program of Higher Education of P. R. China (No.20130043110001), and the Research Fund for Chair Professors of Changbai Mountain Scholars Program.

---