Effective control measures considering spatial heterogeneity to mitigate the 2016–2017 avian influenza epidemic in the Republic of Korea

During the winter of 2016-2017, an epidemic of highly pathogenic avian influenza (HPAI) led to high mortality in poultry and put a serious burden on the poultry industry of the Republic of Korea. Effective control measures considering spatial heterogeneity to mitigate the HPAI epidemic is still a challenging issue. Here we develop a spatial-temporal compartmental model that incorporates the culling rate as a function of the reported farms and farm density in each town. The epidemiological and geographical data of two species, chickens and ducks, from the farms in the sixteen towns in Eumseong-gun and Jincheon-gun are used to find the best-fitted parameters of the metapopulation model. The best culling radius to maximize the final size of the susceptible farms and minimize the total number of culled farms is calculated from the model. The local reproductive number using the next generation method is calculated as an indicator of virus transmission in a given area. Simulation results indicate that this parameter is strongly influenced not only by epidemiological factors such as transmissibility and/or susceptibility of poultry species but also by geographical and demographical factors such as the distribution of poultry farms (or density) and connectivity (or distance) between farms. Based on this result, we suggest the best culling radius with respect to the local reproductive number in a targeted area.


Introduction
During the 2016-2017 winter season, the epidemic of highly pathogenic avian influenza (HPAI) in the Republic of Korea led to high mortality rates in domestic poultry and put a serious economic burden on the poultry industry. By April 4, 2017, 383 farms were reported to be infected by the HPAI virus (subtype H5N6 and H5N8), and approximately 3.7 million poultry (3154 thousand chickens and 332 thousand ducks) from 946 farms were culled (i.e., depopulated) [1]. As most of the culled chickens are layer chickens (2518 thousand), it disrupted the egg supply and led to a surge in the egg price [2,3]. PLOS  Avian influenza (AI) virus is usually not serious in wild birds, but it causes a critical infection in domestic poultry such as chickens and ducks. By ability of the virus to cause disease and mortality in chickens, the infections divided into two forms [4]: low pathogenic avian influenza (LPAI), causing mild symptoms such as decreasing egg production, and the HPAI, causing infected hosts to have high mortality rates up to 90-100% within 2 days [5]. AI is spread by direct contact with birds. The AI virus is mainly transmitted through the feces or respiratory tract of infected birds [6]. Poultry can be infected with the AI virus via betweenfarm movements of vehicles, people, feed and poultry [2,7].
Culling is the process of killing poultry, both infected and uninfected, in areas around the infected region to rapidly contain the spread of an infectious disease with the aim of finally eradicating the disease [8]. Based on a zoning strategy, in practice, authorities mainly carry out following two culling strategies.
• Infected premises (IP) culling: killing all poultry in a farm or small area where the sick or dead poultry are diagnosed as AI-positive.
• Preemptive (PE) culling: preemptively killing all poultry within a protect area, which is declared a area around the IP up to 3 km in the Republic of Korea.
Although such stamping-out policies have been considered to be one of the main control measures for decades, there remain some controversies about animal welfare and effectiveness. Especially, when PE culling is conducted in dense area of farms, it inevitably leads to a huge number of killed poultry. Over 10 million chickens were culled in the Republic of Korea to mitigate H5N1 in 2008 [1]. The decreased poultry industry revenue in China from February to April 2004 caused by HPAI was estimated to 49.6% [9]. Nonetheless, a culling policy is widely adopted because of a dramatic reduction in the risk of dispersion in a short time by decreasing the infectious period of infected farms and also removing susceptible farms which can be infectious in advance.
Mathematical modeling has been used to understand the dynamics of AI epidemics and to provide insights on control measures against the spread of the AI virus. Compartmental SIRtype models without spatial heterogeneity were studied for human infection of the AI virus [10][11][12] with optimal control theory [13], AI transmission between wild birds and poultry [14], and impact of culling strategies on H5N1 infection in domestic bird populations [15,16]. Many studies have explored the role of spatial heterogeneity and control measures in AI epidemics [17][18][19][20][21][22][23]. Agent-based models are beneficial to show warnings on mass culling strategy [17,18]. Distance-dependent transmission probability between poultry farms provides transmission risk in a geographical region [19,24]. Stochastic farm-to-farm transmission models have been used to investigate control measures dependent on spatial heterogeneity, including poultry farm density [19,20]. The transmission characteristics of the HPAI virus were quantified by the estimation of the reproductive number and used it to evaluate the effectiveness of various control measures against HPAI epidemics in recent studies [19,[21][22][23]25].
Even though many mathematical models have focused on various control measures against AI virus, there has been little research on various culling strategies considering the breeding type of poultry farms and the farm density in targeted areas. Moreover, it is still a challenge to suggest an effective control scenario based on the geographical distribution of poultry farms. In this study, we focus on the effects of various culling scenarios on the spread of avian influenza between domestic poultry farms in the Republic of Korea, taking into consideration their breeding types and geographical distributions. We introduce culling functions incorporating the farm density of targeted areas and describe the farm-based metapopulation model with these culling functions. The model is parameterized with AI data reported to the Ministry of Agriculture, Food and Rural Affairs (MAFRA) [1] for the 2016-2017 period. We introduce key parameters to quantify the local transmissibility and use it to compare culling scenarios to maximize the number of surviving farms after the AI epidemic.

Epidemiological and geographical data
In this study, we focus on the 2016-2017 HPAI outbreak in the neighboring sixteen towns which are lower-level administrative divisions of Eumseong-gun and Jincheon-gun in Chungcheongbuk-do. Before the HPAI outbreak, there were 589 poultry farms (454 for chicken and 135 for duck) in these sixteen towns [26]. With the onset of the HPAI outbreak, starting from the town Maengdong, which is located at the center of Eumseong-gun and Jincheon-gun, 71 farms (14 for chicken and 57 for duck) were reported as infected and approximately 2 millions poultry were culled [1] in these sixteen towns.
The data in Table 1 shows the geographical and epidemiological information of the sixteen towns, and is used in our mathematical model. The sequential outbreak data for each type of poultry farms (S1 and S2 Tables) is used for estimation of model parameters. Note that each town is assigned an index in the ascending order of the distance from Maengdong, which is set to index 0. Matrix of pairwise distances between towns (S3 Table) is used in a transmission kernel of our model. Assuming that the farms are uniformly distributed in each town, the farm density is the number of farms per area. The distribution of farm density in a map is shown in Fig 1. The reported numbers of outbreak farms were 57 for duck and 14 for chicken. Maengdong had the largest reported number, i.e., 23, and Iwol-myeon had the second largest number, i.e., 11. The ratio of reported number of farms to total number of farms was 0.12. The ratio of reported duck farms to total number of duck farms was 0.42 while the ratio of reported chicken farms to total number of chicken farms was 0.03. These imply that the 2016-2017 HPAI outbreak occurred mostly at duck farms. Effective control measures considering spatial heterogeneity to mitigate the avian influenza epidemic

Culling rate functions
In this subsection, we describe nonlinear culling rates as a function of the number of AI-confirmed farms. Constant culling rate is used for simplicity in compartmental models [14,16], but it might be inappropriate to consider culling capacity or resource limitation [15,22]. In practice, when a poultry farm is diagnosed positive for HPAI on either clinical diagnosis or laboratory analysis, PE culling is conducted on farms within a protect area which is a declared area around the reported farms of perimeter 3 km in general. The PE culling process is analogous to transmission of diseases; once an infected farm is identified and confirmed AI infection by a epidemiological surveillance, most farms having close contacts would be designated as dangerous contacts and depopulated soon. Therefore, it is reasonable to adopt nonlinear culling terms (or constant culling rates) [15] similar to nonlinear transmission terms such as density-dependent transmission [28,29]. In this work, we employ the nonlinear culling terms in depopulation of the poultry farms. Under the limited culling strategy, we set the culling rate as a rational function of the number of reported farms. Two nonlinear culling types are considered in this study: PE culling and IP culling. PE culling rate function. We set the PE culling rate as a rational function of the number of reported farms [15]. For culling by a reported farm (R), the number of susceptible farm (S) to be culled is proportional to density of those farms within a culling area of radius, r c , i.e., where D is the fraction of culling area, D ¼ pr 2 c =a and a is the area. Then, rate of culling with a saturation factor is given as ZD SR AþR , where η and A are PE culling constant and decay constant for the PE culling rate, respectively. In this study, we define the nonlinear PE culling rate function as cðRÞ ¼ ZD R AþR . Note that ψ(0) = 0, which means PE culling strategy is employed after an AI outbreak occurs.
IP culling rate function. The IP culling process is a control strategy that removes the reported farms to stop disease spreading. Although emergency action guidelines for AI describe that culling procedure should be proceeded in a day [30], it is likely to take 2 days or more because of the simultaneously surging reported cases and limitation of the resources for control measures. We assume that the IP culling rate is a decreasing function of R, and given by �ðRÞ ¼ gB

BþR
, where γ is maximum culling rate and B is a decay constant. When γ is fixed as 1 [30], parameter B can be estimated from the IP culling rates with respect to the number of reported farms by the following approximation: ignoring the type of farm and patch, let N IP be the number of IP culling obtained by data (S4 Table), i.e., where the left term is the cumulative number of AI-reported farms and the right term is the cumulative number of IP culling. These data are provided in S4 Table. In Fig 2, the stars and the solid curve

Metapopulation model with two species
Now we consider a compartmental model to incorporate spatial effects. Let S i , I i , R i denote respectively the number of susceptible, infective but not identified, and reported farms in patch i for i = 0, . . ., 15. Note that in our model I i and R i are infectious, but R i has less transmissibility because of control measures after the case is identified, such as isolation and restriction. Ignoring the species, the spatial-temporal model can be written for i = 0, . . ., 15 as follows The parameter β is the transmission rate, α is the AI virus progression rate to infective farms, and � is the reduction factor of the virus transmissibility for the reported farms. Susceptible farms in town i, S i , can be infected by two ways: either random contacts with infectious farms in the same town, I i and R i , or nonrandom contacts with infectious farms from other towns, I j and R j , for j = 0, . . ., 15 and j 6 ¼ i. Here, contacts between domestic poultry farms represent indirect ways likely to transmit the AI virus via vehicles related to domestic poultry industry. Assuming density-dependent transmission, we use an exponential decay function with respect to the distance between two towns i and j, d ij , as a transmission kernel for nonrandom contact with infectious farms outside the local town i, Kði; jÞ ¼ e À dði;jÞ r 0 , r o is the scaling constant of the transmission kernel. The distance matrix whose the entry in the i-th row and j-th column is d ij is estimated by the shortest path (road) in a map between two towns instead of the Euclidean distance between them. This estimation is reasonable in the case when neighboring two towns have no direct path between them because of some geographical reasons, such as rivers or mountains.
Culling is a localized control measure; it is only implemented in the premises affected by the AI virus. To add this spatial heterogeneity to the model we extend the PE culling rate function. We assume that the susceptible and infective poultry farms in town i are preemptively culled at a rate ψ(R i , D i ; R) when there are AI-reported farms in town i. It is supposed that saturation of culling rate is dependent on the total number of reported farms across all towns, i.e., The PE culling rate, ψ(R i , D i ; R), linearly increases with respect to the number of reported cases in town i once the outbreak occurs, i.e., ψ � ηD i R i , and the fraction of culling area in town i, D i . Finally, the PE culling rate does not continue to increase indefinitely when there exist reported farms in excess. Similarly, we assume that the IP culling is delayed by the total number of reported farms in all towns. Then, we let �ðRÞ ¼ gB BþR with B > 0. Hence the IP culling rate decreases as the total number of reported cases across all towns increases.
We now introduce a metapopulation model with two-type poultry farms: chicken and duck. For i = 0, . . ., 15, S i , I i and R i are divided into two poultry farms such as S ci , I ci , R ci , S di , I di and R di , respectively, where the subscript c is for chicken and d for duck. Then the model with two different farms can be written for i = 0. . ., 15 as the following nonlinear differential equations: where Kði; jÞ ¼ e À dði;jÞ r 0 : The parameter b k 1 k 2 , for k 1 , k 2 2 {c, d}, denotes the transmission rate from k 1 farm to k 2 farm. For example, β cd represents the transmission rate from duck farms to chicken farms. For model simplicity, we suppose that interactions between chicken and duck farms are symmetric, i.e., β cd = β dc . The parameters α c and α d denote the progression rate of chicken and duck farms, respectively. The parameters β cc , β cd , β dd , r 0 , η and B were estimated using the least squares fitting method. The model prediction CðtÞ ¼ R t t 0 cðRðtÞÞIðtÞ þ �ðRðtÞÞRðtÞdt was fitted to the cumulative number of AI-reported farms which is the summation of the cumulative number of PE culled farms that are confirmed as AI-positive and the cumulative number of IP culled farms. The built-in routine lsqcurvefit in MATLAB was used to solve the nonlinear leastsquares problem.
As the first case was reported at a duck farm in Maengdong-myeon during the 2016-17 HPAI epidemic, we set R di (0) = 1 for i = 0 (indicating the infection source town), and R di (0) = 0 for i = 1, . . ., 15. For chicken, there were no initial cases, thereby setting R ci (0) = 0 for all i. Although it was found that there were a few farms already infected when the first case was reported [30], it might be hard to measure the exact number of those farms even by epidemiological surveillance. Therefore, the initial number of infective duck farm I d (0) in Maengdongmyeon was also estimated through the least-squares fitting method.

Reproductive numbers
When an infected individual invades the susceptible population, the average number of secondary infection generated by the primary case over the infectious period, called the basic reproductive number and denoted by R 0 , is an important threshold quantity [31][32][33]. In this work, to find the basic reproductive number, we use the next generation method [33,34]. Let G be the next generation matrix, then R 0 ¼ rðGÞ where ρ is the spectral radius. The (i, j) element of G means how many new infections are introduced into compartment i by the infected from compartment j. We now define the local reproductive number in town j, R ðjÞ 0 , as how many poultry farms are newly infected by infected poultry farms from town j, and it is obtained by the maximum value among the farming types in each town after the sum of each column of G. As we consider two types of poultry farms, the next generation matrix can be written as a 2×2 block matrix, where the block G k 1 k 2 for k 1 , k 2 2 {c, d} is a 16×16 matrix, and the entry of G k 1 k 2 is given by A detailed report can be found in [33].

Parameter estimations
In our study, the progression rate of chicken and duck farms, the initial IP culling rate, and the culling radius were set as α c = 1/2, α d = 1/4, γ = 1, and r c = 3, respectively [1]. We assumed that saturation factor for the PE culling rate, which is the farm number where the culling rate reaches half of its maximum, A = 1 and the infectivity reduction factor, � = 0.01. The model parameters are listed in Table 2.  myeon (Index 0). The fitted model agrees well with the observed data for both types of farms. Note that the transmission rate between duck farms (0.00707, 95% confidence interval: 0.00484-0.00930) is approximately 1.6 times higher than the transmission rate between chicken farms (0.00441, 95% confidence interval: 0.00183-0.00700). Interestingly, even though we only estimated non-spatial parameters such as transmission rates and scaling constants in the metapopulation model with two species, the total number of AI-positive farms from the model demonstrated a good fit to the corresponding data in most towns; the root-mean-square errors between data and model for chicken, and duck and total farms of the sixteen towns are 0.9392, 1.5020 and 1.5572, respectively.

Reproductive numbers
Based on the initial population of poultry farms in Table 1 and the estimated parameters in Table 2, the estimate R 0 for the AI outbreak was 1.3427.    (2)) with all other parameters in Table 2. The left panel depicts the final size of susceptible farms with respect to culling radius. The right panel displays the total number of PE culling (dot-dashed), IP culling (dashed) and both (solid) with respect to culling radius. A larger culling radius increases the loss by PE culling and suppresses outbreaks thereby decreases the loss by IP culling. This implies that the two culling strategies are in a compensatory relationship. The total loss of farms by culling reaches its minimum at r c best = 2.24 km of culling radius in which the final size of the AI epidemic reaches its peak.

PE culling strategy
Although we found the best culling radius, r c best , at which the total loss of farms by culling reaches its minimum, spatial heterogeneity was not considered. We set culling radius as localized parameter; r c high , culling radius of towns with R ðiÞ 0 > 1 and r c low , culling radius of towns with R ðiÞ 0 < 1. Fig 6 depicts the impact of the two localized culling radii on the final size of susceptible farms by a color bar. In the parameter domain, the final size of susceptible farms along the axis r c high dramatically increases at first, then slowly decreases later. Meanwhile, the final size of susceptible farms decreases as r c low increases. The final size of susceptible farms is maximized as 459 when r c high and r c low are 2.65 and 0, respectively. The best culling radius in towns with R ðiÞ 0 > 1, r c high = 2.65 km, is smaller than the government's policy (3 km), but greater than one ignoring the spatial heterogeneity, r c best = 2.24 km.

Discussion
Using farm-to-farm transmission dynamics incorporating the two poultry types, chicken and duck, we estimated the spread of the HPAI outbreak in the Republic of Korea in 2016-2017.  We found that from modeling result the transmissibility between duck farms was higher than that between either chicken farms or different type of farms. Ducks can carry and shed the AI virus without symptoms and have low mortality while chickens have high pathogenicity and mortality [5]. In addition, the outbreak started at a duck farm in high farm density area. These epidemiological and spatial factors amplified the virus and allowed it to spread easily surrounding poultry farms. Therefore, duck farms played an important role in the spread of the HPAI virus in Eumseong-gun and Jincheon-gun. This is the first study on estimating the basic reproductive number of the 2016-2017 HPAI epidemic in the Republic of Korea and introducing the local reproductive number as an indicator of the virus transmission in a given area. There were six towns in which the local reproductive number is greater than 1: Maengdong (R ð0Þ 0 ¼ 1:6095), Daeso (R ð3Þ 0 ¼ 1:5144), Deoksan (R ð1Þ 0 ¼ 1:3913), Iwol (R ð8Þ 0 ¼ 1:2940), Samseong (R ð10Þ 0 ¼ 1:0926), and Jincheon (R ð7Þ 0 ¼ 1:0277). These towns seem to be either relatively close to Maengdong (Deoksan, Geumwang and Daeso) or have high farm density (Maengdong, Deoksan and Samseong). The distance from Maengdong to Deoksan, Geumwang and Daeso are 7.76 km, 9.77 km and 9.95 km, respectively, which are almost half of mean distance (18.31 km) to Maengdong. Since the transmission rate depends on the distance between towns by the kernel, Kði; jÞ ¼ e À dði;jÞ r 0 , it is clear that the local reproductive number is affected by the distance. The high farm density in Deoksan and Samseong (1.31 and 1.60, respectively) might allow the virus to spread easily surrounding poultry farms via movement of humans (farm personnel and visitors) and vectors (rodents), or contaminated environment (air and water) [35].
Maengdong, which is the town with the highest local reproductive number, is the region with the highest number of duck farms and the source of outbreak. The three towns with the Effective control measures considering spatial heterogeneity to mitigate the avian influenza epidemic lowest local reproductive number are Gamgok, Soi, and Eumseong, which also have the lowest farm density. These imply that farm density plays an important role on the 2016-2017 HPAI epidemic in the Republic of Korea, which is in line with the previous works on livestock diseases [19,20,[36][37][38][39][40]. Furthermore, we found that R ðiÞ 0 is strongly influenced by not only epidemiological factors, such as the number of duck farms with high transmissibility in town i, but also by geographical factors, such as the location of town i and its proximity to other towns with high susceptibility.
Our finding showed that culling strategies need to take heterogeneity into account for reducing loss of poultry farms. We suggested that PE culling has to be focused on the area in which the local reproductive number is greater than 1, and the culling radius must be greater than that, ignoring the spatial heterogeneity; Reinforced PE culling in the area wherein the local reproductive number is greater than 1 could allow for early depopulation of infected farms before infection spreads to dense area. Meanwhile, in the areas in which the local reproductive number is less than 1, losses from the infection might be greater than losses from PE culling.

Conclusions
In this paper, we presented mathematical modeling for the spatial-temporal transmission dynamics of HPAI using geographical and epidemiological data of the 2016-2017 AI epidemic in the Republic of Korea. To consider spatial heterogeneity, we introduced the PE culling rate as a function of the number of reported farms [15] with different coefficients based on the farm density in each town. As the exact location of the poultry farms affected by the AI virus was not available, a metapopulation framework has been adopted with the two types of poultry farms assuming random mixing between farms in a patch and nonrandom contact by the transmission kernel [19,41] between farms in different patches. We introduced the local reproductive numbers using the next generation method [33,34] in the metapopulation modeling framework. This quantified parameter allowed one to assess the transmissibility in each town. For total and even both types of farms, the model predictions and data about the AI-reported farms are in good agreements. The estimated transmission rates showed that the transmission between duck farms played an important role in the spread of the HPAI virus in Eumseong-gun and Jincheon-gun.
Our findings showed that the local reproductive number could be an indicator of the likelihood of virus transmission in a given area. It was revealed that this parameter was strongly influenced not only by epidemiological factors such as transmissibility and/or susceptibility of poultry species but also by geographical and demographical factors such as the distributions of poultry farms (or density) and connectivity (or distance) between farms. Based on this result, we found that the culling radius of PE culling should be adjusted by considering the local reproductive number in the target area. Therefore, to determine which area is supposed to be more strictly controlled during an AI outbreak, we believe that veterinary/public health officials can use the local reproductive number in a real-time warning system. Our study can be applied to other animal diseases (e.g., foot-and-mouth disease [42], brucellosis [25]) in which such heterogeneities are crucial factors to be considered.
Supporting information S1