Transmission dynamics of brucellosis: Mathematical modelling and applications in China.

Brucellosis, the most common zoonotic disease worldwide, represents a great threat to animal husbandry with the potential to cause enormous economic losses. Meanwhile, brucellosis is one of the major public-health problems in China, and the number of human brucellosis cases has increased dramatically in recent years. In order to show the main features of brucellosis transmission in China, we give a systematic review on the transmission dynamics of brucellosis including a series of mathematical models and their applications in China. For different situations, dynamical models of brucellosis transmission in single population and multiple populations are devised based on ordinary differential equations. Furthermore, we revealed the spatial-temporal characteristics and effective control measures of brucellosis transmission. The results may provide new perspectives for the prevention and control of other types of zoonoses.

Introduction

Brucellosis, a bacterial disease caused by various Brucella species, is one of the most common zoonotic infections in the world [1], [2], [3]. There are four Brucella species which are mainly responsible for the transmission of the disease: B. melitensis in goats and sheep, B. abortus typically found in cattle, B. canis in dogs, and B. suis in swine [4], see Fig. 1. Even though these four species of Brucella can infect humans, B. melitensis remains the major cause of human disease worldwide (and may account for up to 90% of all brucellosis cases [5]). The remaining illnesses are caused by B. abortus and B. suis, with rare but persisting cases of B. canis infections in humans [6]. There are three main transmission ways. People can be infected through eating undercooked meat and unpasteurized dairy products both which are carrying brucella. Meanwhile, people get infected by inhaling brucella, mainly laboratory workers who work with the bacteria. In addition, can also infect workers who are in close contact with animals or animal waste can be infected by the skin, wounds or mucous membranes, containing slaughterhouse workers, meat-packing plant employees, veterinarians and even hunters. Although there are a small number of reports of vertical and horizontal transmission between humans [7], [8], it is generally acknowledged that human-to–human transmission of the infection is a very rare event [9].

Fig. 1

Main animal hosts of brucella and different main human hosts, where the main susceptible hosts is in the light blue box. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

According to the length and severity of symptoms, the disease in humans is arbitrarily classified as acute (less than 8 weeks), subacute (from 8 to 52 weeks), or chronic (more than 1 year) [9]. The disease is commonly underestimated, misdiagnosed and once chronic disease develops resistant to treatment, which consists of antibiotics for long periods [10], [11], [12]. Mortality is reported to be negligible, but the illness can persist for several years. Though less severe in animals, brucellosis can cause great economic losses by adversely affecting reproduction and fertility, survival of newborns, and milk yields [13], [14].

Human brucellosis matches the regions with high levels of animal infection endemicity in the world: the Mediterranean basin, Middle East, Western Asia, Africa, and South America [3] where hundreds of thousands of new cases are reported annually [9], [15]. The developed countries have been successfully controlled the animal brucellosis. However, many of these control options are less achievable in developing countries [16]. In China, human brucellosis is a class B notifiable infectious disease, and each confirmed case must be reported to the Chinese CDC (CCDC) through the National Notifiable Disease Surveillance System (NNDSS) since 2004 [17]. Although many measures based on the control programs for brucellosis have been set up, the brucellosis-positive rate in humans has increased significantly in recent years [18], [19]. The annual and cumulative reported human brucellosis cases of mainland China and the selected eleven provinces with high incidence rate from 2004 to 2018 are listed in Table 1. The spatiotemporal distribution map presented in Fig. 2 showed that human brucellosis was widely distributed in the northern, northeastern, and western China. On the contrary, cases are more sporadic in other places. At present, brucellosis epidemic areas in China can be divided into three categories: (1) Severe epidemic areas: Heilongjiang, Jilin, Inner Mongolia, Ningxia, Qinghai, Xizang, Xinjiang, Henan; (2) General epidemic areas: Shaanxi, Gansu, Sichuan, Hebei, Liaoning, Shandong, Shanxi, Guangdong, Guangxi; (3) Distributed epidemic areas: other provinces and cities. Meanwhile, the spatial spread of brucellosis in China is mainly from pastoral areas to rural and urban areas in recent years, the prevalence rate among livestock in some areas has exceeded that in pastoral areas. On the other hand, the spread of human Brucellosis in space shows a gradual spread from east to west and from north to south over time. Therefore, it is of great significance to study the transmission mechanism and characteristics of brucellosis in different epidemic areas and effective control measures, which has become a public health urgent need to solve the problem.

Table 1

Reported cases of human brucellosis in mainland China and the 11 provinces with highest rates during 2004–2018.

Region	2004	2005	2006	2007	2008	2009	2010	2011
Mainland China	11472	18416	19013	19721	27767	35816	33772	38151
Inner Mongolia	4356	8663	7951	8117	11105	16551	16224	17320
Shanxi	1599	2320	3452	4040	4834	4768	3888	5135
Heilongjiang	2219	3943	2949	2377	3670	4724	4861	5178
Hebei	718	1181	2334	2404	3173	3218	2503	3149
Xinjiang	317	431	357	386	484	525	914	1477
Jilin	444	521	599	877	2139	3452	2905	2063
Henan	234	141	194	295	593	703	791	1225
Liaoning	1039	604	439	343	484	509	606	853
Shaanxi	265	483	636	741	1034	911	525	556
Shandong	41	163	108	138	135	156	222	449
Ningxia	1	26	7	28	94	157	207	471
Region	2012	2013	2014	2015	2016	2017	2018
Mainland China	39515	43486	57222	56989	47139	39026	21217
Inner Mongolia	12017	8911	10142	7238	5970	7451	5227
Shanxi	6130	6895	8586	6997	4587	3088	1485
Heilongjiang	7459	7056	5628	5960	5342	4298	2292
Hebei	4066	5103	6526	5526	3774	2733	1621
Xinjiang	2476	4468	7567	8820	8402	6076	2869
Jilin	1956	2024	1811	1625	1480	1206	568
Henan	1926	3256	5200	5573	3993	2526	1031
Liaoning	1475	2053	2811	2922	2338	1927	1198
Shaanxi	596	795	1475	1221	948	730	368
Shandong	665	1278	2809	3691	3886	3198	1411
Ningxia	449	883	2064	2888	2160	1660	1059

Fig. 2

Spatiotemporal distribution of annual human brucellosis cases, by province of China in 2004–2018.

Mathematical models have the potential to analyze the mechanisms of transmission and complex epidemiological characteristics of infectious diseases, and can put forward new approaches to prevent and control future epidemics. In 2010, Ainseba et al. [20] proposed two ways of brucellosis transmission in sheep: direct contact with infected individuals and indirect exposure to the contaminated environment. By mathematical analysis, they obtained the net reproduction number of dynamical model. Based on the sheep population in Algeria, they numerically simulated the different cases of direct infection rate and indirect infection rate and got the impact of two transmission modes on the disease. In addition, they provided an annual slaughtering policy and affirmed the importance of indirectly spread in the persistence of the disease. In recent years, our team has focused on the study of brucellosis transmission in different endemic areas of brucellosis in mainland China through the establishment of dynamics models and proposed corresponding control measures. In severe epidemic areas, Li et al. [21] proposed a deterministic model in Hinggan League of Inner Mongolia to describe the transmission of brucellosis for sheep. We estimated the controlled reproduction number to be about 1.9789 and concluded that vaccination, detection and the mixed cross infection between basic ewes and other sheep are the very important factors for brucellosis transmission in Hinggan League. Based on the mathematical analysis, we provided the effective measures, which is combination of prohibiting mixed feeding between basic ewes and other sheep, vaccination, detection and elimination. Besides, Nie et al. [22] researched the impact of transmission about dairy cattle brucellosis considering the direct mode and indirect mode in Jilin province. The results indicated importing, culling and sterilizing were very important factors for brucellosis in Jilin province. In distributed epidemic areas, Zhang et al. [23] discussed the brucellosis transmission in Zhejiang province and concluded environment transmission was major factor to arouse the prevalence of the epidemic and the epidemic would be a periodic phenomenon in Zhejiang province. According to the dynamical analysis, self-supplying production of the dairy cows and sterilization times of twice a week in the infected regions are the most effective measures to control brucellosis transmission.

This paper is organized as follows. In the next section, we show the present situation and transmission dynamics of brucellosis combined different mathematical dynamic models. In Section 3, spatial–temporal characteristics analysis and simulation results about brucellosis will be given. In addition, we will propose some disease control strategies. Section 4 demonstrates a brief discussion about main results, shortcomings and future work.

Transmission dynamical models of brucellosis

The dynamical models play an critical role in epidemic forecast and control, not only for anthroponosis and zoonosis, but also for anthropozoonosis and vector borne diseases. The greatest advantage of dynamical models is based on propagation mechanism of epidemic to study its propagation tendency under different assumptions and precondition, with the limitation of data material and lack of experimental conditions. Therefore, it can be seen that dynamical models are the most important method and play a significant role in any epidemics forecasting and strategy making. Dorigatti et al. applied a dynamical model with a spatial transmission kernel to analyze the between-farm transmission of the H7N1 highly pathogenic avian influenza virus [26]. For Severe acute respiratory syndrome (SARS) outbroken in 2003, Lipsitch applied dynamical method to estimate the reproduction number and infectious cases for SARS [27]. Hao et al. established a dynamical model to present the full-spectrum dynamics of COVID-19 in Wuhan [28]. Mukandavire et al. applied dynamical model to study transmission dynamics of cholera and assess the magnitude of necessary interventions to control epidemic disease [29]. In the review of foot and mouth disease, Keeling introduced three different models which were used to predict the disease dynamics and provide control measures for its outbreak in the UK [30].

The dynamical method is to establish different types of differential equation or difference equation to describe the variation of dependent variables with independent variables. The independent variable is general time, space or age. Three key elements during the course of epidemic development are susceptible groups, source of transmission and route of transmission. Then, susceptible host population, infectious host population, vector and so on can be considered as dependent variables, respectively. The route of transmission can be exhibited as the interaction between dependent variables.

Basic dynamical model

Taking the total host population for example, whose number at time t is denoted by , its variation with time can be given as the following literal equation.

The change rate of the total host population at time t = its increment rate at time t - its decrement rate at time t.

The change rate of the total host population at time t can be expressed as derivative of with respect to time, denoted by . The increment rate at time t is denoted by , and the decrement rate at time t is denoted by . The schematic diagram of change of can be seen in Fig. 3. Then we can transform above literal equation into the following differential equation:

Fig. 3

The flow chart of model (2.1).

According to the developing process of brucellosis, the host population can be divided into the most basic subpopulations: the susceptible host population (also called as healthy host population), the infected host population, whose numbers with time t can be denoted by and , respectively. Assume that the increment population of the host population is all healthy, then the increment rate at time t only occurs in susceptible group and it is assumed to be constant A. The decrement rate is generally assumed to be linear related to the total population, , where d is the decrement rate coefficient. The propagation process of brucellosis among host population can be thought as the transmission of brucella caused by random contact between individuals, which can be described as a transmission rate in the model. Assume that the contact number per individual per unit time is , where is the number of infected individuals among contact group. Then the contact number of all susceptible individuals with infectious individuals per unit time is . The infection probability after one contact is denoted by , then the number of new infected individuals per unit time is . When is linear related to the total population . Denoting , bilinear incidence ratio is obtained, which is used for the case that the contact number per individual per untie time is proportional to the total number of individuals [21], [24], [31], [32], [22], [34], [36], [37], [40], [41]. When is constant, . Denoting , standard incidence ratio is obtained, which is used for the case that the contact number per individual per untie time is constant [23], [35].

According to above assumptions, the schematic diagram of the transmission of brucella between subpopulations S and I can be seen in Fig. 4 and corresponding simple SI brucellosis model with bilinear incidence rate is given as equations:

Fig. 4

The flow chart of model (2.2).

Dynamical models of brucellosis transmission in single population

Bacteria is present in host animals’ fluids (urine, milk, placental fluid, etc.), and is a significant cause of abortions and sterility in these animals. Therefore, the infection can occur among hosts through direct contact with infected animals or their contaminated animal products, or indirectly by contact with animals’ fluids and contaminated water, soil, food, air and dust in the environment. Multiple transmission routes are considered in researches [21], [24], [25], [31], [23], [32], [22], [33], [34], [35], [36], [37], [40], [41], [42]. The transmission of environment-to-individual in most works is described by the form of bilinear incidence, , where is the transmission rate coefficient and W is the amount of bacteria in environment. The saturating incidence rate is adopted in [34], which describes the saturation effects of environment transmission when the amount of bacteria in environment increases to some degree, and this kind of transmission term is more reasonable than bilinear or standard incidence rate and is closer to actual situation. Considering the existence of minimal infecting dose invading into animal individuals, Zhang et al. gave a switch incidence ratio to describe indirect transmission of environment to individuals [42]. Li et al. proposed a more general incidence rate [25], [35]. In addition, for the brucella in environment, its variation includes the discharge of brucella into environment by infected individuals I, generally denoted by and is the discharge load per individual per unit time, and its decay, denoted by lW and l is decay rate coefficient.

Infection of brucellosis may be acute, chronic (causing symptoms for years), or may be asymptomatic. The incubation period of brucellosis is 1 to 3 weeks, up to several months, with an average of 2 weeks. Animals with latent infection do not show any clinical symptoms or signs, but play an important role in the spread of the epidemic. Most researches consider the exposed population as a separate compartment, generally denoted by E, to establish the dynamical model. The incidence ratio is described as , where is the transmission ratio coefficient of the infectious subpopulation I to the susceptible subpopulation S and is the transmission ratio coefficient of the exposed subpopulation E to the susceptible subpopulation S.

With the influence of climate on transmissibility and survival of bacteria, and culture pattern of domestic animals, the outbreak of brucellosis has obvious seasonality and the periodic propagation mechanism which is described in references [33], [43]. In the previous work, we took month as unit of time and adopted periodic parameters in model to describe periodic infectivity, the birth of animal, the sale of animal and abortion time [33]. Beauvais et al. adopted switching model in accordance with season to investigate the spread of brucellosis on a mixed sheep-and-cattle farm [43].

Combining above characteristics, we give the more general flow Chart of propagation of brucella in a single population, which is divided into three subpopulations: and I, and their survival environment W, see Fig. 5. The corresponding nonautonomous model is given in the following system:where and are periodic infection rate coefficient, is periodic sale rate coefficient, is quantity of bacterium discharged by each animal at time is the clinical outbreak rate coefficient, is culling rate coefficient of the infectious population I due to disease and l is the decay rate coefficient of brucella in environment.

Fig. 5

The flow chart of model (2.3). Black solid line represents flow of the individuals and decay of brucella in environment, black dotted line represents the discharge of brucella into environment, and blue dotted line represents the infection of brucella in environment to the susceptible. The variation of susceptible population S includes the birth or import of individuals which are assumed to be healthy, the sale of susceptible individuals, and the being infected of susceptible individuals to enter the exposed population. The variation of exposed population E includes entering of susceptible individuals being infected, the sale of exposed individuals, and showing clinical symptoms of exposed individuals to enter the infected population. The variation of infected population I includes entering of exposed individuals showing clinical symptoms, the sale of infected individuals and the culling of infected individuals. The variation of the brucella in environment W includes the discharge of brucella into environment by exposed individuals and infected individuals, and its decay. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Brucella are predominantly transmitted via direct contact with abortion and birth fluids and semen of infected animals, so the adult animals are mainly infection resource. Moreover, lambs who are less than 6 months are less likely to be infected, and adult fattened sheep are easily infected. The age structure is considered in [36], [37], [40], [43]. Li and Hou proposed a multi-stage dynamic model, involving young sheep population and adult sheep population, to describe sheep brucellosis transmission [37], [40]. Beauvais et al. considered age categories as follows: young which cannot become infected, juvenile which can become infected but not infectious, and adult which can become infectious when they have a late abortion or give birth [43].

Dynamical models of brucellosis transmission in multiple populations

Infected animals become infectious and carry bacteria to spread among the same species in the common pasture, farm or vehicle, and then can be transmitted to humans through skin and mucous membrane contact, digestive tract, respiratory tract and other channels. Human-to–human transmission is very rare. Then the transmission among animals is crossed and bidirectional, which is described by multiple population dynamical models [21], [24], [25], [44], [45], [46]. The cross-infection between i species and j species is described by and , where represents the number of the susceptible in i species, represents the number of the infected in j species, represents the transmission rate coefficient of to . Other parameters have similar meanings, not given in detailed. We give the dynamical model of brucellosis transmission between 2 species in the following equations [24] (the corresponding flow chart in Fig. 6):

Fig. 6

The flow chart of model (2.4) with two animals populations. The subpopulations of the 1st species are denoted by and . The subpopulations of the 2nd species are denoted by and . The brucella in their environment are denoted by and , respectively. Besides the variation of subpopulations in Fig. 5, the variation of susceptible population and adds the vaccination of susceptible individuals to enter the immune subpopulation and the entering of immune individuals since vaccine protection has expired. The variation of immune population and adds the entering of susceptible individuals due to vaccination and being susceptible of immune individuals since vaccine protection has expired. In this model, the susceptible population can be infected not only by the exposed , the infected and the brucella in environment , but also by the exposed , the infected and the brucella in environment . The figure is reproduced from [24].

Each population is divided into four subpopulations: the susceptible S, the exposed E, the infected I, the immune subpopulation R and their survival environment V. is the clinical outbreak rate coefficient. is the cull rate coefficient. is the discharge rate coefficient of brucella into environment. is the decay rate coefficient of brucella in environment. represents the immunization rate coefficient of i species and represents the immune lost rate coefficient of vaccine.

The transmission of animal to human is unidirectional, which kind of transmission mechanism is described in researches [21], [25], [32], [33], [34], [35], [36], [47]. In this case, the spread of the epidemic is only determined by the animal population, so the equations about animals are independent of equations about human and can be analyzed separately:which considered the brucellosis spread between sheep and human and the flow chart can be seen in Fig. 7. Sheep population is divided into two subpopulations: the susceptible and the infected . After infected with brucellosis, the infection stage of human is divided into acute stage and chronic stage. Acute stage refers to the case that appears the clinical symptoms and the course of disease is within 6 months. Chronic stage refers to the case that the disease cannot be cured after 6 months. After standardized treatment, 90% of the patients in the acute stage can be cured. Some human cases can become chronic which is difficult to cure. They will face long-term recurrent attacks, physical weakness, be difficult to work and live normally. So, in model (2.5) [35], human population is divided into four subpopulations: the susceptible , the exposed and the acute infected , and the chronic infected . Their survival environment is denoted by W. is the standard incidence rate of the infected animal to susceptible human. is the incidence rate of the brucella in environment to susceptible human. The meaning of other parameters can be referred to the work in Ref. [35]. In Ref. [48], the transmission of wildlife to domestic animal is considered to be unidirectional.

Fig. 7

The flow chart of model (2.5) between human and sheep populations. In this model, the birth number of sheep population is in proportion to the total sheep population. b is birth rate coefficient, that is the birth number of each sheep per unit time, and is equal to the sale rate coefficient. Similarly, the birth number of human population is in proportion to the total human population. is birth rate coefficient of human population and is death rate coefficient of human population. is the transformation rate coefficient from the exposed to the acute infected. p is the transformation rate coefficient from the acute infected to the chronic infected. The figure is reproduced from [35].

Patch dynamical models of brucellosis transmission

With the development and convenience of transportation, the inter-provincial transport and trade of domestic animals has resulted in propagation of brucellosis over long distance. Then the prevalence of brucellosis in multiple regions is interactional and should be discussed in detailed, and the population in each region is divided into several subclasses according to infection status. The works [49], [50] adopted patch model to describe the infection of the epidemic among discrete regions by applying different methods. The model in [49] considered the immigration of animal among n patches (see Fig. 8) and the immigrant from ith patch into jth patch will belong to jth patch:withwhose flow chart can be seen in Fig. 9. denote the immigration rate of the susceptible, the exposed, the infectious, and brucella in environment from jth patch to ith patch, respectively. denote the emigration rate of the susceptible, the exposed, the infectious, and brucella in environment in ith patch, respectively. They satisfy the following relationships in (2.7). The meaning of other parameters can be referred to [49].

Fig. 8

The schematic diagram of the migration of individuals. Here, the red solid line represents the migration of individuals among patches. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9

The flow chart of model (2.6) with patches. The red lines represent the immigration of individuals and brucella among patches. In i-th patch, the subpopulations are denoted by and , and the brucella in its environment is denoted by . In this model, the spread of the disease among individuals within the patch is the same as Fig. 6 except having no consideration of immune subpopulation. Besides, the infectivity of the exposed to the susceptible is less than the infectivity of the infected, which is multiplied by . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In contrast, the model in [50] considered the cross infection among patches and the subpopulation of ith patch that contacts with the population of jth patch still belongs to ith patch, and the model is similar to multi-population cross-infection model. In addition, Cantrell proposed a reaction diffusion model with the Laplace operator to describe the spread of the epidemic throughout continuous space [51].

Other dynamical models

In China, current prevention and control measures mainly include periodic testing, immunization and quarantine for domestic animals, publicity measures for farmers, sterilization for environment and so on, where testing, immunization and sterilization are assessed by applying dynamical method in [21], [23], [32], [22], [33], [34], [35], [40], [46], [52], [53]. These measures are described as the linear term in model, such as , where k is the immunization rate coefficient, is the detection rate coefficient and l is the decay rate coefficient of brucella in environment. In addition, the involved domestic animals will be tested according to the report of infected farmers. This measure can lead to the checkout and culling of part of positive individuals, which is called as tracking-culling measure in [34] by establishing culling term , where I is the number of infected animals, H is the number of susceptible human, U is the number of the infected human and p is the number of family members in each household.

In fact, the transmission of the epidemic will be affected by all kinds of random factors. We established a stochastic dynamical model to describe the random import of the infected cattle [23]:where are the proportion of the exposed and the infected about the import number A, respectively. is a time series of random deviates derived from the normal distribution with mean zero and unit variance, which is stochastic perturbations to and , and represents the intensity of . The meaning of other parameters can be referred to [23]. Besides, some researchers applied stochastic simulation model [54], [55] to describe transmission mechanism and assess the effect of control measures.

Applications of mathematical modelling in China

Spatial-temporal characteristics analysis

Spatial-temporal statistical methods have been used to study the transmission characteristics of infectious diseases in space and time, and have achieved remarkable results [57], [58], [59], [60]. In the following part, the spatial-temporal characteristics of brucellosis cases in China are analyzed, including spatial distribution characteristics, global spatial autocorrelation, cluster and outlier and hot spot analysis.

Spatial distribution characteristics

We collect the brucellosis cases data to draw geographic visualization map and analyze its spatial epidemiological characteristics. Fig. 10 shows the spatial distribution of cumulative brucellosis cases in every city of China from 2010 to 2018. The brucellosis cases have spread across all the provinces in mainland China, mainly distributed in the North, Northwest and Northeast China, including most cities of Inner Mongolia autonomous Region, Heilongjiang, Jilin, Liaoning, Hebei, Shanxi, Shaanxi, Shandong, Henan, Ningxia Hui Autonomous Region and Xinjiang Uygur Autonomous Region. By comparing with the annual spatial distribution of brucellosis cases as shown in Fig. 11, we conclude that brucellosis distributed in 28 provinces, 210 cities of China in 2010. By 2013, brucellosis have spread to all provinces in China, even in the cities of South and Southeast China. The brucellosis cases in the cities in Inner Mongolia autonomous Region and its surrounding provinces are obviously severer than in other areas during this period. And the number of cases in these areas began to decline since 2012.

Fig. 10

Geographic distribution of cumulative brucellosis at city level in China from 2010 to 2018. Colored gradients reflect the number of brucellosis cases. The brucellosis cases have spread across all the provinces in mainland China, mainly in the North, Northwest and Northeast China.

Fig. 11

Geographic distribution of brucellosis cases at city level in China per year from 2010 to 2018. Colored gradients reflect the number of cases. The brucellosis cases in Inner Mongolia autonomous Region and its surrounding provinces are obviously higher than in other areas during this period. And the number of cases in these areas began to decline since 2012.

Global spatial autocorrelation

Spatial autocorrelation analysis is a spatial statistical method to test the attribute value of an element as to whether it is dramatically associated with its adjacent unit [61]. Spatial autocorrelation includes global and local spatial autocorrelation [62]. Global spatial autocorrelation is used to measure the overall clustering tendency in the study region, evaluate whether the pattern expressed is clustered, dispersed, or random. In this study, we use global Moran’s I statistic to investigate spatial association regarding brucellosis cases at city level in China. The range of Moran’s I value is . At the level of statistical significance, Moran’s implies positive autocorrelation, and a higher positive Moran’s I implies that the distribution is more aggregated; Moran’s indicates negative autocorrelation, and a lower negative Moran’s I indicates that the distribution is more dispersed; and Moran’s means no autocorrelation, suggesting that the disease is randomly distributed [63]. The result is listed in the Table 2, it shows that the Z-score is greater than 1.96, the P-values are significant at 0.05 every year, which demonstrates that the brucellosis cases are not randomly distributed and that there is global autocorrelation among brucellosis cases at city level in China. And the Global Moran’s I value is greater than zero, which indicates that it is positive autocorrelation and the distribution of brucellosis cases is spatially clustered.

Table 2

The spatial autocorrelation test of brucellosis cases at city level in China, 2010–2018.

Year	Moran’s I	Expected I	Variance	Z-score	P-value	Result
2010	0.121054	-0.002755	0.002227	2.623678	0.008699	clustered
2011	0.136076	-0.002755	0.002268	2.915200	0.003555	clustered
2012	0.162996	-0.002755	0.002297	3.458230	0.000544	clustered
2013	0.192833	-0.002755	0.002330	4.052223	0.000051	clustered
2014	0.200052	-0.002755	0.002341	4.191201	0.000028	clustered
2015	0.225382	-0.002755	0.002418	4.775494	0.000002	clustered
2016	0.232085	-0.002755	0.002404	4.652856	0.000003	clustered
2017	0.198830	-0.002755	0.002397	4.117383	0.000038	clustered
2018	0.185647	-0.002755	0.002209	4.008537	0.000061	clustered

Cluster and outlier analysis

Local Moran’s I can be further used to clarify the patterns and the exact location of the clusters among local areas [64]. At the level of statistical significance, Moran’s implies a spatial clustering of similar values surrounding the unit area i; Moran’s indicates a spatial clustering of non-similar values surrounding the unit area i. A statistically significant cluster area will be represented with four patterns of spatial correlation:high-high clusters (HH); low-low clusters (LL); high-low clusters (HL); and low–high clusters (LH). In reality, the high-high clusters (i.e.high infection areas surrounded by other high infection areas) are the most important patterns for disease prevention and control [65]. The Local Moran’s I statistics is performed to describe the distribution patterns of brucellosis cases at city level in China, which is shown in the Fig. 12, Fig. 13.

Fig. 12

Local Indicators of Spatial Association (LISA) maps for cumulative brucellosis cases at city level in China from 2010 to 2018; H-H, high-high cluster; L-L, low-low cluster; L-H, low–high cluster; H-L, high-low cluster. There is only one distribution pattern of the brucellosis cases, namely high-high cluster (HH) (green color), which are mainly centered in the most cities of Inner Mongolia autonomous Region, a few cities of its surrounding provinces and Xinjiang Uygur Autonomous Region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 13

Local Indicators of Spatial Association (LISA) maps for brucellosis cases at city level in China per year, from 2010 to 2018; H-H, high-high cluster; L-L, low-low cluster; L-H, low-high cluster; H-L, high-low cluster. There are two distribution patterns of the brucellosis cases, namely high-high cluster (HH) (green color) and low-high (LH) outlier (blue color). The HH cluster occurs every year, is mainly concentrated in most cities of Inner Mongolia Autonomous Region and a few cities of its surrounding provinces. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 shows the cluster and outlier locations of cumulative brucellosis cases at city level in China from 2010 to 2018. There is only one distribution pattern of the brucellosis cases, namely high-high cluster (HH) (green color), which are mainly centered in most cities of Inner Mongolia autonomous Region; a few cities of Heilongjiang, Jilin, Hebei, Shanxi and Shaanxi provinces around Inner Mongolia Autonomous Region; and Ili Kazak Autonomous Prefecture, Karamay City and Altay Region of Xinjiang Uygur Autonomous Region. It indicates that the number of brucellosis cases is relative high in these cities and their surrounding cities. Fig. 13 shows the cluster and outlier locations of brucellosis cases at city level in China per year from 2010 to 2018. There are two distribution patterns of the brucellosis cases, namely high-high cluster (HH) (green color) and low-high (LH) outlier (blue color). The HH cluster occurs every year, and are mainly concentrated in most cities of Inner Mongolia Autonomous Region and a few cities of surrounding provinces. It also has been distributed in several cities of Xinjiang Uygur Autonomous Region since 2013, and the HH cluster area in 2015 and 2016 are wider. The LH outlier (blue color) distribution only occurs in one or two cities of Xinjiang Uygur Autonomous Region in 2015, 2016 and 2017, which means the number of brucellosis cases in the city is lower than that of the surrounding cities. The number of brucellosis cases of other cities is zero or shows no statistically significant data and random distributions.

Hot spot analysis

Getis-Ord can be used to recognize the distribution of hot spot areas and cold spot areas, that is, to determine the high infection areas and low infection areas [66]. The Getis-Ord statistic itself is the statistical test of Z-score. At the level of statistical significance, implies hot spot areas, and a higher positive Z-score implies that the clustering of high-value is more intense; indicates cold spot areas, and a lower negative Z-score indicates that the clustering of low-value is more loosened; and Z-score closes to 0 means no spatial clustering, suggesting that the disease is randomly distributed. We use Getis-Ord statistics to describe the hot spot distribution of brucellosis cases at city level in China, which is shown in the Fig. 14, Fig. 15.

Fig. 14

Hot spot analysis using Getis-Ord for cumulative brucellosis cases at city level in China, from 2010 to 2018. There are two concentrated hot spot areas (). One hot spot area mainly located in most cities of Inner Mongolia autonomous Region and a few cities of its surrounding provinces, the other one located in some cities of Xinjiang Uygur Autonomous Region.

Fig. 15

Hot spot analysis using Getis-Ord for brucellosis cases at city level in China per year, from 2010 to 2018. The cities of Inner Mongolia autonomous Region and its surrounding provinces occurred hot spots about the epidemic of brucellosis every year from 2010 to 2018. Since 2014, some cities of Xinjiang Uygur Autonomous Region started to appear hot spots pattern.

Fig. 14 shows the hot and cold spot areas of cumulative brucellosis cases at city level in China from 2010 to 2018. There are two concentrated hot spot areas (), which indicated these areas have a relatively high number of brucellosis cases that formed a spatial high value cluster distribution. One hot spot area mainly located in most cities of Inner Mongolia autonomous Region, and some cities of Heilongjiang, Jilin, Liaoning, Hebei and Shanxi, all of which are around Inner Mongolia Autonomous Region. The other concentrated hot spot areas are located in some cities of Xinjiang Uygur Autonomous Region. Fig. 15 shows the hot and cold spot areas of brucellosis cases at city level in China per year, from 2010 to 2018. The cities of Inner Mongolia autonomous Region and its surrounding provinces appeared hot spots about the epidemic of brucellosis every year from 2010 to 2018. Since 2014, some cities in Xinjiang Uygur Autonomous Region started to appear hot spots pattern. It means that these areas have a relatively high number of brucellosis cases which formed a spatial cluster distribution. These hot spot areas and their surrounding areas need to strengthen prevention and control.

Transmission dynamics and control measures of brucellosis in China

Mathematical modeling has the potential to analyze the mechanisms of transmission and the complexity of epidemiological characteristics of infectious diseases, and can provide new approaches to prevent and control future epidemics [56]. In recent years, several mathematical dynamic modeling studies have reported the transmission of brucellosis with actual infected cases in China [35], [21], [23].

Human brucellosis in China

Human brucellosis matches the regions of the world with high levels of animal infection endemicity: the Mediterranean basin, Middle East, Western Asia, Africa, and South America where hundreds of thousands of new cases are reported annually [35]. In mainland China, sheep and goats probably were the main animal hosts transmitting the disease to humans [24]. Although many control measures of brucellosis have been set up, the brucellosis-positive rate in humans has increased significantly in recent years. Li et al. devised model (2.5) to assess the fitting data of the cumulative brucellosis cases of mainland China and eleven provinces with high case numbers in Table 1 with Monte Carlo simulation, and the plots are shown in Fig. 16 [35]. The parameter values of model (2.5) are showing the following section.

Fig. 16

Brucellosis model fitting for the cumulative incidence of human brucellosis cases in mainland China and eleven selected provinces. The dotted red line represents the model simulation, light grey shaded area shows the 95% confident interval (CI) for all 1000 simulations, and the blue circles mark the reported data for cumulative human brucellosis cases. The figure is reproduced from [35]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We know the period of human brucellosis is about two weeks, so the clinical outcome rate for exposed people is per year. From the China Statistical Yearbook [38] and the China Animal Husbandry Statistical Yearbook [39], one can obtain the demographic parameter values (including human populations birth rate and death rates , sheep recruitment and slaughter rate b) for mainland China and the 11 provinces with the highest incidence for human brucellosis listed in Table 3. Li et al. [35] fix the human indirect transmission rate , and assume that for mainland China and 11 selected provinces with the highest incidence for human brucellosis. By using extensive Markov-chain Monte-Carlo (MCMC) simulations based on the adaptive combination Delayed rejection and Adaptive Metropolis (DRAM) algorithm, the estimates of and for these 11 provinces and mainland China can be obtained, and are shown in Table 4.

Table 3

Values of and and 95% confidence intervals ().

	b 95% CI	b_h 95% CI (10^-3)	d_h 95% CI (10^-3)
Mainland China	0.9026 (0.8531–0.9521)	12.098 (11.9863–12.2097)	6.937 (6.7422–7.1318)
Xingjiang	0.8363 (0.7470–0.9256)	15.918 (15.5573–16.2787)	4.973 (4.7345–5.2125)
Shandong	1.3203 (1.2255–1.4151)	11.682 (11.3864–11.8776)	6.353 (6.1953–6.5107)
Liaoning	0.9914 (0.9723–1.0104)	6.382 (6.0943–6.6697)	5.758 (5.3973–6.1187)
Henan	1.0794 (1.0524–1.1064)	11.616 (11.4145–11.8162)	6.459 (6.3734–6.6166)
Ningxia	0.8515 (0.7890–0.9139)	14.509 (13.7663–15.2517)	4.755 (4.5832–4.9268)
Shanxi	0.5087 (0.4892–0.5283)	11.201 (10.7591–11.6429)	5.819 (5.6489–5.9891)
Hebei	1.2397 (1.0895–1.3899)	12.91 (12.649–13.171)	6.544 (6.2962–6.6918)
Heilongjiang	0.7864 (0.7605–0.8123)	7.448 (7.1854–7.7106)	5.538 (5.2718–5.8042)
Shannxi	0.6855 (0.6541–0.717)	10.115 (9.9321–10.2979)	6.1516 (6.0809–6.2321)
Inner Mongolia	1.0228 (1.0021–1.0435)	9.525 (9.262–9.7878)	5.7025 (5.5515–5.8535)
Jilin	0.7667 (0.7292–0.7956)	7.2026 (6.6134–7.7994)	5.2854 (5.0902–5.4805)

Table 4

Estimated values of and with their 95% confidence intervals ().

	β_s 95% CI	β_sw 95% CI	β_h 95% CI
Mainland China	0.5583 (0.5522–0.5643)	0.1125 (0.1086–0.1186)	0.0676 (0.0662–0.0691)
Xinjiang	0.2395 (0.1602–0.3188)	0.3037 (0.2832–0.3243)	0.7291 (0.6525–0.8057)
Shandong	0.5739 (0.3944–0.7534)	0.3743 (0.3319–0.4176)	0.9223 (0.7376–1.1070)
Liaoning	0.1870 (0.1138–0.2602)	0.2970 (0.2816–0.3123)	0.0738 (0.0297–0.1180)
Henan	0.2389 (0.1829–0.2950)	0.3402 (0.3115–0.3689)	0.4457 (0.3881–0.5032)
Ningxia	0.3761 (0.2102–0.5420)	0.2847 (0.2455–0.3240)	0.4264 (0.2148–0.6380)
Shanxi	0.3878 (0.3622–0.4134)	0.0544 (0.0437–0.0650)	0.9355 (0.9049–0.9661)
Hebei	0.6213 (0.5381–0.7045)	0.1847 (0.1550–0.2143)	0.9295 (0.8963–0.9627)
Heilongjiang	0.3781 (0.3176–0.4386)	0.1311 (0.1167–0.1456)	0.5256 (0.4561–0.5951)
Shaanxi	0.4598 (0.4023–0.5173)	0.0579 (0.0427–0.0731)	0.0929 (0.0331–0.1565)
Inner Mongolia	0.6275 (0.6030–0.6520)	0.1437 (0.1344–0.1530)	0.4638 (0.4401–0.4874)
Jilin	0.5347 (0.3661–0.7034)	0.1984 (0.1657–0.2311)	0.4628 (0.3405–0.5851)

The basic reproduction number provides useful guidelines for the prevention and control strategies of epidemics. The basic reproduction number of model (2.5) iswhere and are partial reproduction numbers due to environment-to-individual transmission and individual-to-individual transmission, respectively. By using parameter values of model (2.5), they obtained that the basic reproduction numbers of eleven selected provinces and mainland China (see Table 5). Table 5 shows that the local basic reproduction numbers of provinces with an obvious increase in incidence were much larger than the average of the whole country.

Table 5

Estimated values of and and 95% confidence intervals.

	R0i	95% CI	R0e	95% CI	R_0	95% CI
Mainland China	0.6191	(0.6043–0.6338)	0.5198	(0.4566–0.5826)	1.1389	(1.0608–1.2164)
Xinjiang	0.2863	(0.2102–0.3623)	1.5136	(1.3940–1.6327)	1.7999	(1.6042–1.9950)
Shandong	0.5622	(0.4266–0.6979)	1.0326	(0.8294–1.2361)	1.5948	(1.2560–1.9340)
Liaoning	0.2908	(0.1804–0.4011)	1.1356	(0.8990–1.3722)	1.4264	(1.0793–1.7734)
Henan	0.2473	(0.2008–0.2939)	1.2847	(1.1658–1.4039)	1.5319	(1.3665–1.6978)
Ningxia	0.5473	(0.1861–0.9084)	1.2718	(0.8519–1.6911)	1.8190	(1.0381–2.5995)
Shanxi	0.6741	(0.5805–0.7674)	0.5365	(0.4349–0.6372)	1.2106	(1.0154–1.4047)
Hebei	0.6398	(0.5999–0.6869)	0.4766	(0.4440–0.5089)	1.1163	(1.0439–1.1957)
Heilongjiang	0.7906	(0.6960–0.8852)	0.3741	(0.2904–0.4583)	1.1646	(0.9863–1.3435)
Shaanxi	0.6874	(0.5964–0.7784)	0.3671	(0.2839–0.4504)	1.0545	(0.8802–1.2288)
Inner Mongolia	0.6817	(0.6332–0.7302)	0.5145	(0.4220–0.6074)	1.1962	(1.0552–1.3376)
Jilin	0.7616	(0.4946–1.0287)	0.8271	(0.6326–1.0222)	1.5887	(1.1272–2.0509)

Animals vaccination (v), quarantine, separation and elimination of the infected animals for brucellosis (), and the disinfection of environment (l is the disinfection frequency) can be considered as control measures for brucellosis. Li et al. concluded the minimum vaccination coverage rate, removal rate, and disinfection frequency which are required to control brucellosis epidemics for the eleven selected provinces and mainland China in Table 6 [35]. They suggested that a combination of animal vaccination, environment disinfection, and elimination of infected animals would be necessary to ensure cost-effective control for brucellosis.

Table 6

Estimates of minimum vaccination coverage rate v, removal rate , and disinfection frequency l.

	Vaccination rate	Removal rate	Disinfection frequency
Mainland China	0.1487	0.1254	3
Xinjiang	0.5420	0.6689	9
Shandong	0.4548	0.7853	10
Liaoning	0.3646	0.4277	5
Henan	0.4235	0.5742	6
Ningxia	0.5491	0.6974	14
Shanxi	0.2121	0.1071	5
Hebei	0.1271	0.1442	3
Heilongjiang	0.1724	0.1295	6
Shaanxi	0.0630	0.0374	2
Inner Mongolia	0.2000	0.2007	5
Jilin	0.4519	0.4514	18

In previous paper [21], Li et al. proposed a deterministic model to describe the transmission of brucellosis among sheep and from sheep to humans in Hinggan League of Inner Mongolia. They used the deterministic model to fit the annually infected human brucellosis data in Hinggan League from 2001 to 2011 (see Fig. 17 (a)), which indicated that their model provided a good match to the reported data. By analyzing the tendency of newly infected human brucellosis cases in next 14 years (see Fig. 17 (b)), newly infected human brucellosis cases would increase in the next seven or eight years, and then decrease slowly. They used 1000 samples of Latin hypercube sampling method to sample the derived expression for control reproduction number of Hinggan League, the distribution in values is shown for in Fig. 18. From the method of Latin hypercube sampling about these parameters, they obtained that the mean value of is 1.9789, the minimum value of is found to be 1.2648 and the maximum value of is 5.039.

Fig. 17

(a) The comparison between the reported human brucellosis data in Hinggan League of Inner Mongolia from 2001 to 2011 and the simulation of newly infected human brucellosis cases from the model. (b) The tendency of newly infected human brucellosis data in Hinggan League of Inner Mongolia in next 14 years since 2011. The figure is reproduced from [21].

Fig. 18

Histogram of the estimated aggregate value from Latin hypercube sampling. The figure is reproduced from [21].

In Hinggan League, basic ewes and other sheeps are often mixed feeding together. Therefore, there exist single species internal infection from basic ewes or other sheep, and the mixed cross infection for brucellosis between basic ewes and other sheeps. In order to investigate the influence of prohibiting mixed feeding between basic ewes and other sheep, Li et al. considered the variations of newly infected human brucellosis cases for different control parameters with mixed feeding and without mixed feeding (see Fig. 19) [21]. Comparing two figures of Fig. 19, they concluded when the efficient vaccination rate and the seropositive detection rate of basic ewes and other sheep reach , the control reproduction number is about , the final scale of newly infected human brucellosis cases of Hinggan League will tend to zero in the future. But in Fig. 19 a, when the efficient vaccination rate and the seropositive detection rate of basic ewes and other sheep reach , the control reproduction number is about , which means Hinggan League brucellosis will become endemic. Hence, the mixed cross infection between basic ewes and other sheep plays an important role in the persistence of brucellosis, it is a very important factor for brucellosis transmission. They concluded that combination of prohibiting mixed feeding between basic ewes and other sheep, vaccination, detection and elimination were the useful macro-control strategies in controlling human brucellosis in Hinggan League.

Fig. 19

The variations of newly infected human brucellosis cases for different control parameters (vaccination rate , detection rate c and elimination rate ) since 2011. (a) With mixed feeding, , , respectively. (b) Without mixed feeding, , , respectively. The figure is reproduced from [21].

Animal brucellosis in China

Brucellosis is a notifiable disease that can infect cows, swine, goats, sheep, dogs and humans. Brucellosis infection of animal can bring about serious economic consequences for the animal husbandry. Not only does it cause losses from abortion, but it also creates marketing limitations. When an infection case occurs, the disease can quickly spread by the herd, causing significant productivity losses. In the Zhejiang province, brucellosis between dairy cows has attracted significant attention of the public and government, which is a major public-health and economically devastating zoonosis. Hence, Zhang et al. proposed an SEIV dynamical model to investigate the internal transmission dynamics of brucellosis in dairy herd, predicted the infection situation, and assessed the prevention and control measures [23]. The prediction trends of brucellosis in Zhejiang province by deterministic model and stochastic model are shown in Fig. 20. With random import of infected dairy cows, the prediction situations can have large differences which can reach up to several hundreds, and the disease in dairy cows cannot be controlled without taking more effective measures. Due to constraint of sensibility of detection methods, as long as there are imported dairy cows, the input of infected dairy cows is unavoidable. Under certain circumstances, the epidemic would become a periodic phenomenon [23].

Fig. 20

100 times of prediction results of with deterministic model and stochastic model during 50 years. The figure is reproduced from [23].

Discussion and conclusion

Human beings and animals are supposed to coexist peacefully in the same environment, which is named as “One Health” [67]. However, due to international trade of exotic animals, and the increasingly frequent human encroachment on wildlife habitats, coupled with the development of international travel networks and urbanization, have disrupted the boundaries of human animal environment. As we known, zoonoses not only affect people health, but also have influences on the structure of the whole society. Brucellosis is the most common zoonosis in China and thus we give a systematic review on the transmission dynamics and control measures of brucellosis in China in this review. We constructed a series of dynamical models for different situations including human brucellosis and animal brucellosis and data fitting is used to obtain the appropriate parameters value. Moreover, quantitative assessment of control measures containing vaccination, detection and elimination for different regions in China are presented based on these results.

It should be pointed out that this review is mainly based on deterministic models of brucellosis transmission. However, the mechanisms of transmission may be affected by precipitation amount, temperature change, and other climatic factors [68], [69], [70]. Under these circumstances, stochastic models are the better methods to portray the propagation mechanism of the spread of brucellosis. Meanwhile, it should also be noted that spatial effects play an important role in brucellosis transmission, which may form pattern formation with rich structures [71], [72], [73]. For such reason, one may need to combine the method of transmission dynamics based on partial differential equations (mainly reaction–diffusion equations) [74] and geographic information system (GIS) to associate the outbreak and development trend of brucellosis with the specific geographical location, and build a brucellosis prediction and early warning system. Additionally, big data analysis is also required to reveal the features of brucellosis transmission in both space and time [75], [76]. The ultimate goal is to build a framework for key technologies for prediction and early warning signal of brucellosis transmission. These issues will be systematically studied in the future investigations.

We want to pose the future topics on brucellosis transmission in three directions. Firstly, cluster degree is an important feature to describe the spatial evolution of infectious diseases [77], [78]. Consequently, how to combine cluster degree with dynamical models for describing the spread of brucellosis is a fundamental problem to be addressed. Secondly, economy is closely related to zoonoses and thus one needs to find out which types of control measures are optimal based on the economy cost and efficiency of disease control [79]. Lastly, transport is one of the most important links for Livestock trading and thus how to quantify the effects of transport between different cities or countries is an urgent question to be solved [80].

As we known, COVID-19 is one of the hottest events which can not be avoided when you look at what happened in the last year [81], [82]. As of August 13th, 2020, there are about 20 million people infected by COVID-19 and 700 thousand people died all over the world due to it. COVID-19 is also one of the zoonosis which warns people not to excavate nature endlessly. In order to prevent the emergence and spread of zoonosis, and reduce its harm to society, wildlife experts, veterinary experts, sociologists, political scientists, anthropologists, economists and other relevant experts should work together with interdisciplinary and international cooperation.

CRediT authorship contribution statement

Gui-Quan Sun: Data curation, Formal analysis, Investigation, Methodology, Methodology, Validation, Writing - original draft, Writing - review & editing. Ming-Tao Li: Data curation, Formal analysis, Investigation, Writing - original draft, Writing - review & editing. Juan Zhang: Data curation, Formal analysis, Investigation, Writing - original draft, Writing - review & editing. Wei Zhang: Investigation, Writing - original draft. Xin Pei: Data curation, Writing - original draft. Zhen Jin: Methodology, Validation, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.