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

Jonggul Lee; Youngsuk Ko; Eunok Jung

doi:10.1371/journal.pone.0218202

Abstract

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.

Citation: Lee J, Ko Y, Jung E (2019) Effective control measures considering spatial heterogeneity to mitigate the 2016–2017 avian influenza epidemic in the Republic of Korea. PLoS ONE 14(6): e0218202. https://doi.org/10.1371/journal.pone.0218202

Editor: Daniel Becker, Montana State University, UNITED STATES

Received: February 13, 2019; Accepted: May 28, 2019; Published: June 13, 2019

Copyright: © 2019 Lee et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The data that support the findings of this study are publicly available online at www.mafra.go.kr, www.data.go.kr and www.kosis.kr. The authors confirm that the data set used to reach the conclusions drawn in the manuscript are available in the manuscript and Supporting Information files.

Funding: This work is resulted from the Konkuk University research support program. This paper is supported by the Korea National Research Foundation (NRF) grant funded by the Korean government (MEST) (NRF-2017R1A2B2004651). This paper is also supported by National Institute for Mathematical Sciences (NIMS) grant funded by the Korea government (MSIT) (No. B19610000).

Competing interests: The authors have declared that no competing interests exist.

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].

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 between-farm 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 SIR-type models without spatial heterogeneity were studied for human infection of the AI virus [10–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–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–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.

Materials and methods

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.

Download:

Fig 1. Map of Eumseong-gun and Jincheon-gun, which are the municipal-level divisions located in the central part of the Republic of Korea, including the density of poultry farms with white indicating low density and red indicating high density.

The index of the town is written below its name. There are 16 towns; 9 towns in Eumseong-gun and 7 in Jincheon-gun. Reprinted from [27] under a CC BY license, with permission from National Geographic Information Institute, original copyright 2017.

https://doi.org/10.1371/journal.pone.0218202.g001

Download:

Table 1. Numbers of poultry farms and reported cases of sixteen towns in 2016-2017.

https://doi.org/10.1371/journal.pone.0218202.t001

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, and a is the area. Then, rate of culling with a saturation factor is given as , 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 . 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 , 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., N_IP = ϕ(R)R, then ϕ(R) = N_IP/R. By integrating dR/dt from (2), we have , 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 denote the approximated IP culling rate with respect to the number of AI-reported farms and the corresponding calibrated curve, respectively. Fig 2 shows that the IP culling rate decreases as the number of reported farms increases. Note that the average of IP culling rate from the data might be overestimated when there are a large number of AI-reported farms.

Download:

Fig 2. Estimation of IP culling rate.

The stars and the solid curve denote the approximated IP culling rate with respect to the number of AI-reported farms and the corresponding calibrated curve, respectively. The dashed line represents constant culling rate averaged over data. The data for IP culling rate approximation is presented in S4 Table.

https://doi.org/10.1371/journal.pone.0218202.g002

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 (1) 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 ≠ 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, , 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., . Then, we let with A > 0. 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_iR_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 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: (2) The parameter , for k₁, k₂ ∈ {c, d}, denotes the transmission rate from k₁ farm to k₂ 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₀, η and B were estimated using the least squares fitting method. The model prediction 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 least-squares 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 Maengdong-myeon 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 , is an important threshold quantity [31–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 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, , 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, (3) where the block for k₁, k₂ ∈ {c, d} is a 16×16 matrix, and the entry of is given by (4) A detailed report can be found in [33].

Results

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.

Download:

Table 2. Model parameters.

https://doi.org/10.1371/journal.pone.0218202.t002

Fig 3 displays the cumulative number of AI-positive farm (circles) and the best-fitted model results (solid curves) after the first case was reported on November 16, 2016, in Maengdong-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.

Download:

Fig 3. In each panel, circles represent the AI-reported farm data and black solid curves represent the corresponding simulation results.

https://doi.org/10.1371/journal.pone.0218202.g003

Reproductive numbers

Based on the initial population of poultry farms in Table 1 and the estimated parameters in Table 2, the estimate for the AI outbreak was 1.3427. Fig 4 displays the local reproductive number of each town by a map and graph which shows the correlation among local reproductive number, the farm density and the number of duck farms. In Fig 4, the towns with the local reproductive number greater than 1 are Maengdong (), Daeso (), Deoksan (), Iwol (), Samseong (), and Jincheon (). In most towns, the local reproductive numbers are positively correlated to the farm density (with correlation coefficient, ρ = 0.6606). Furthermore, the local reproductive numbers are also positively correlated to the number of duck farms (with correlation coefficient, ρ = 0.9383).

Download:

Fig 4.

(A) Local reproductive numbers in each town. The darker red color represents the larger reproductive number. (B) Local reproductive number with respect to farm density in each town. The size of circle shows the size of duck farms. Reprinted from [27] under a CC BY license, with permission from National Geographic Information Institute, original copyright 2017.

https://doi.org/10.1371/journal.pone.0218202.g004

PE culling strategy

Fig 5 shows the impact of culling radius on the AI outbreak using the metapopulation model (Eq (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.

Download:

Fig 5. Impact of culling radius on the AI outbreak using the metapopulation model (Eq (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. Note that when culling radius is 2.24 km, the final size of susceptible farms is maximized and the total number of culling is minimized.

https://doi.org/10.1371/journal.pone.0218202.g005

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 and r_{c low}, culling radius of towns with . 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_{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.

Download:

Fig 6. Impact of local-dependent culling radii on the final size of susceptible farms.

The final size of susceptible farms with respect to culling radii is maximized as 459 when r_{c high} and r_{c low} are 2.65 and 0, respectively.

https://doi.org/10.1371/journal.pone.0218202.g006

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 (), Daeso (), Deoksan (), Iwol (), Samseong (), and Jincheon (). 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, , 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 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–40]. Furthermore, we found that 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 Table. Time series outbreak data (chicken farms) of the sixteen towns [1].

https://doi.org/10.1371/journal.pone.0218202.s001

(CSV)

S2 Table. Time series outbreak data (duck farms) of the sixteen towns [1].

https://doi.org/10.1371/journal.pone.0218202.s002

(CSV)

S3 Table. Pairwise distance matrix of the sixteen towns (km).

https://doi.org/10.1371/journal.pone.0218202.s003

(CSV)

S4 Table. Daily AI report and IP culling number.

This are gathered and translated from HPAI infection daily reports of the Ministry of Agriculture, Food and Rural Affairs [1].

https://doi.org/10.1371/journal.pone.0218202.s004

(CSV)

References

1. Ministry of Agriculture, Food and Rural Affairs. HPAI situation report;. Available from: http://www.mafra.go.kr.
2. Animal and Plant Quarantine Agency. Animal disease control;. Available from: http://www.qia.go.kr.
3. The United States Department of Agriculture. Korea—Republic of, Poultry and Products Annual Report; 2017.
4. Centers for Disease Control and Prevention, National Center for Immunization and Respiratory Diseases (NCIRD). Avian Influenza in Birds;. Available from: https://www.cdc.gov/flu/avianflu/avian-in-birds.htm.
5. Kang Y, Shen X, Yuan R, Xiang B, Ren T. Pathogenicity and transmissibility of three avian influenza A (H5N6) viruses isolated from wild birds. Journal of Infection. 2018;76(3):286–294. pmid:29307740
- View Article
- PubMed/NCBI
- Google Scholar
6. Vandegrift KJ, Sokolow SH, Daszak P, Kilpatrick AM. Ecology of avian influenza viruses in a changing world. Annals of the New York Academy of Sciences;1195(1):113–128. pmid:20536820
- View Article
- PubMed/NCBI
- Google Scholar
7. jin Lee H, Suh K, su Jung N, bok Lee I, hwan Seo I, kyung Moon O, et al. Prediction of the spread of highly pathogenic avianinfluenza using a multifactor network: Part2- Comprehensive network analysis with direct/indirect infection route. Biosystems Engineering. 2014;118(1):115–127.
- View Article
- Google Scholar
8. World Organisation For Animal Health. Avian Influenza Portal;. Available from: http://www.oie.int/en/animal-health-in-the-world/avian-influenza-portal/prevention-and-control/.
9. Elci C. The Impact of HPAI of the H5N1 Strain on Economies of Affected Countries. Izmir University of Economics; 2006. p. 104–117. Available from: https://ideas.repec.org/h/izm/prcdng/200610.html.
10. Iwami S, Takeuchi Y, Liu X. Avian-human influenza epidemic model. Mathematical Biosciences. 2007;207(1):1–25. pmid:17010999
- View Article
- PubMed/NCBI
- Google Scholar
11. Iwami S, Takeuchi Y, Korobeinikov A, Liu X. Prevention of avian influenza epidemic: What policy should we choose? Journal of theoretical biology. 2008;252:732–741. pmid:18384816
- View Article
- PubMed/NCBI
- Google Scholar
12. Iwami S, Takeuchi Y, Liu X. Avian flu pandemic: Can we prevent it? Journal of theoretical biology. 2009;257:181–190. pmid:19094996
- View Article
- PubMed/NCBI
- Google Scholar
13. Jung E, Iwami S, Takeuchi Y, Jo Tc. Optimal control strategy for prevention of avian influenza pandemic. Journal of Theoretical Biology. 2009;260(2):220–229. pmid:19501107
- View Article
- PubMed/NCBI
- Google Scholar
14. Pandit PS, Bunn DA, Pande SA, Aly SS. Modeling highly pathogenic avian influenza transmission in wild birds and poultry in West Bengal, India. Scientific Reports. 2013;3(1):2175. pmid:23846233
- View Article
- PubMed/NCBI
- Google Scholar
15. Gulbudak H, Martcheva M. Forward hysteresis and backward bifurcation caused by culling in an avian influenza model. Mathematical Biosciences. 2013;246(1):202–212. pmid:24056242
- View Article
- PubMed/NCBI
- Google Scholar
16. Shim E, Galvani AP. Evolutionary repercussions of avian culling on host resistance and influenza virulence. PLoS ONE. 2009;4(5):1–8.
- View Article
- Google Scholar
17. Kim T, Hwang W, Zhang A, Sen S, Ramanathan M. Multi-agent model analysis of the containment strategy for Avian Influenza (AI) in South Korea. Proceedings—IEEE International Conference on Bioinformatics and Biomedicine, BIBM 2008. 2008; p. 335–338.
18. Kim T, Hwang W, Zhang A, Sen S, Ramanathan M. Multi-agent modeling of the South Korean avian influenza epidemic. BMC Infectious Diseases. 2009;10(5):236.
- View Article
- Google Scholar
19. Boender GJ, Hagenaars TJ, Bouma A, Nodelijk G, Elbers ARW, De Jong MCM, et al. Risk maps for the spread of highly pathogenic avian influenza in poultry. PLoS Computational Biology. 2007;3(4):704–712.
- View Article
- Google Scholar
20. Backer JA, van Roermund HJW, Fischer EAJ, van Asseldonk MAPM, Bergevoet RHM. Controlling highly pathogenic avian influenza outbreaks: An epidemiological and economic model analysis. Preventive Veterinary Medicine. 2015;121(1-2):142–150. pmid:26087887
- View Article
- PubMed/NCBI
- Google Scholar
21. Bett B, Henning J, Abdu P, Okike I, Poole J, Young J, et al. Transmission rate and reproductive number of the H5N1 highly pathogenic avian influenza virus during the december 2005-July 2008 epidemic in Nigeria. Transboundary and Emerging Diseases. 2014;61(1):60–68. pmid:22925404
- View Article
- PubMed/NCBI
- Google Scholar
22. Stegeman A, Bouma A, Elbers A, de Jong M, Nodelijk G, de Klerk F, et al. Avian Influenza A Virus (H7N7) Epidemic in The Netherlands in 2003: Course of the Epidemic and Effectiveness of Control Measures. The Journal of Infectious Diseases. 2004;190(12):2088–2095. pmid:15551206
- View Article
- PubMed/NCBI
- Google Scholar
23. Ward MP, Maftei D, Apostu C, Suru A. Estimation of the basic reproductive number (R0) for epidemic, highly pathogenic avian influenza subtype H5N1 spread. Epidemiology and Infection. 2009;137(2):219–226. pmid:18559127
- View Article
- PubMed/NCBI
- Google Scholar
24. Dorigatti I RRPABL Mulatti P. Modelling the spatial spread of H7N1 avian influenza virus among poultry farms in Italy. Epidemics. 2010;2(1):29–35. pmid:21352774
- View Article
- PubMed/NCBI
- Google Scholar
25. Darbon A, Valdano E, Poletto C, Giovannini A, Savini L, Candeloro L, et al. Network-based assessment of the vulnerability of Italian regions to bovine brucellosis. Preventive Veterinary Medicine. 2018;158(March):25–34. pmid:30220393
- View Article
- PubMed/NCBI
- Google Scholar
26. Open Data Porta. Eumsung-gun domestic animals statistical research;. Available from: http://www.data.go.kr.
27. National Geographic Information Institute. Administrative boundaries;. Available from: http://data.nsdi.go.kr/dataset/20180927ds0059.
28. Begon M, Bennett M, Bowers RG, French NP, Hazel SM, Turner J. A clarification of transmission terms in host-microparasite models: numbers, densities and areas. Epidemiology and infection. 2002;129(1):147–153. pmid:12211582
- View Article
- PubMed/NCBI
- Google Scholar
29. McCallum H, Barlow N, Hone J. How should pathogen transission be modeled? Trends in Ecology and Evolution. 2001;16(6):295–300.
- View Article
- Google Scholar
30. Ministry of Agriculture, Food and Rural Affairs. HPAI emergency action guidelines;. Available from: http://www.mafra.go.kr.
31. Heffernan JM, Smith RJ, Wahl LM. Perspectives on the basic reproductive ratio. Journal of the Royal Society, Interface / the Royal Society. 2005;2(4):281–293.
- View Article
- Google Scholar
32. van den Driessche P. Reproduction numbers of infectious disease models. Infectious Disease Modelling. 2017;2(3):288–303. pmid:29928743
- View Article
- PubMed/NCBI
- Google Scholar
33. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences. 2002;180(1):29–48. pmid:12387915
- View Article
- PubMed/NCBI
- Google Scholar
34. Diekmann O, Heesterbeek JAP, Metz J. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology. 1990;28(4):365–382. pmid:2117040
- View Article
- PubMed/NCBI
- Google Scholar
35. Assessing the probability of introduction and spread of avian influenza (AI) virus in commercial Australian poultry operations using an expert opinion elicitation. PLoS ONE. 2018;13(3):1–19.
- View Article
- Google Scholar
36. Gilbert M, Pfeiffer DU. Risk factor modelling of the spatio-temporal patterns of highly pathogenic avian influenza (HPAIV) H5N1: A review. Spatial and Spatio-temporal Epidemiology. 2012;3(3):173–183. pmid:22749203
- View Article
- PubMed/NCBI
- Google Scholar
37. Hill SC, Lee YJ, Song BM, Kang HM, Lee EK, Hanna A, et al. Wild waterfowl migration and domestic duck density shape the epidemiology of highly pathogenic H5N8 influenza in the Republic of Korea. Infection, Genetics and Evolution. 2015;34:267–277. pmid:26079277
- View Article
- PubMed/NCBI
- Google Scholar
38. Manore C, McMahon B, Fair J, Hyman JM, Brown M, Labute M. Disease properties, geography, and mitigation strategies in a simulation spread of rinderpest across the United States. Veterinary Research. 2011;42(1):1–12.
- View Article
- Google Scholar
39. Meadows AJ, Mundt CC, Keeling MJ, Tildesley MJ. Disentangling the influence of livestock vs. farm density on livestock disease epidemics. Ecosphere. 2018;9(7).
- View Article
- Google Scholar
40. Tildesley MJ, Bessell PR, Keeling MJ, Woolhouse MEJ. The role of pre-emptive culling in the control of foot-and-mouth disease. Proceedings Biological sciences / The Royal Society. 2009;276(1671):3239–48.
- View Article
- Google Scholar
41. Backer JA, Hagenaars TJ, Nodelijk G, van Roermund HJW. Vaccination against foot-and-mouth disease I: Epidemiological consequences. Preventive Veterinary Medicine. 2012;107(1-2):27–40. pmid:22749763
- View Article
- PubMed/NCBI
- Google Scholar
42. Buhnerkempe MG, Tildesley MJ, Lindström T, Grear DA, Portacci K, Miller RS, et al. The impact of movements and animal density on continental scale cattle disease outbreaks in the United States. PLoS ONE. 2014;9(3):1–10.
- View Article
- Google Scholar

[ref1] 1. Ministry of Agriculture, Food and Rural Affairs. HPAI situation report;. Available from: http://www.mafra.go.kr.

[ref2] 2. Animal and Plant Quarantine Agency. Animal disease control;. Available from: http://www.qia.go.kr.

[ref3] 3. The United States Department of Agriculture. Korea—Republic of, Poultry and Products Annual Report; 2017.

[ref4] 4. Centers for Disease Control and Prevention, National Center for Immunization and Respiratory Diseases (NCIRD). Avian Influenza in Birds;. Available from: https://www.cdc.gov/flu/avianflu/avian-in-birds.htm.

[ref5] 5. Kang Y, Shen X, Yuan R, Xiang B, Ren T. Pathogenicity and transmissibility of three avian influenza A (H5N6) viruses isolated from wild birds. Journal of Infection. 2018;76(3):286–294. pmid:29307740
View Article
PubMed/NCBI
Google Scholar

[6] View Article

[7] PubMed/NCBI

[8] Google Scholar

[ref6] 6. Vandegrift KJ, Sokolow SH, Daszak P, Kilpatrick AM. Ecology of avian influenza viruses in a changing world. Annals of the New York Academy of Sciences;1195(1):113–128. pmid:20536820
View Article
PubMed/NCBI
Google Scholar

[10] View Article

[11] PubMed/NCBI

[12] Google Scholar

[ref7] 7. jin Lee H, Suh K, su Jung N, bok Lee I, hwan Seo I, kyung Moon O, et al. Prediction of the spread of highly pathogenic avianinfluenza using a multifactor network: Part2- Comprehensive network analysis with direct/indirect infection route. Biosystems Engineering. 2014;118(1):115–127.
View Article
Google Scholar

[14] View Article

[15] Google Scholar

[ref8] 8. World Organisation For Animal Health. Avian Influenza Portal;. Available from: http://www.oie.int/en/animal-health-in-the-world/avian-influenza-portal/prevention-and-control/.

[ref9] 9. Elci C. The Impact of HPAI of the H5N1 Strain on Economies of Affected Countries. Izmir University of Economics; 2006. p. 104–117. Available from: https://ideas.repec.org/h/izm/prcdng/200610.html.

[ref10] 10. Iwami S, Takeuchi Y, Liu X. Avian-human influenza epidemic model. Mathematical Biosciences. 2007;207(1):1–25. pmid:17010999
View Article
PubMed/NCBI
Google Scholar

[19] View Article

[20] PubMed/NCBI

[21] Google Scholar

[ref11] 11. Iwami S, Takeuchi Y, Korobeinikov A, Liu X. Prevention of avian influenza epidemic: What policy should we choose? Journal of theoretical biology. 2008;252:732–741. pmid:18384816
View Article
PubMed/NCBI
Google Scholar

[23] View Article

[24] PubMed/NCBI

[25] Google Scholar

[ref12] 12. Iwami S, Takeuchi Y, Liu X. Avian flu pandemic: Can we prevent it? Journal of theoretical biology. 2009;257:181–190. pmid:19094996
View Article
PubMed/NCBI
Google Scholar

[27] View Article

[28] PubMed/NCBI

[29] Google Scholar

[ref13] 13. Jung E, Iwami S, Takeuchi Y, Jo Tc. Optimal control strategy for prevention of avian influenza pandemic. Journal of Theoretical Biology. 2009;260(2):220–229. pmid:19501107
View Article
PubMed/NCBI
Google Scholar

[31] View Article

[32] PubMed/NCBI

[33] Google Scholar

[ref14] 14. Pandit PS, Bunn DA, Pande SA, Aly SS. Modeling highly pathogenic avian influenza transmission in wild birds and poultry in West Bengal, India. Scientific Reports. 2013;3(1):2175. pmid:23846233
View Article
PubMed/NCBI
Google Scholar

[35] View Article

[36] PubMed/NCBI

[37] Google Scholar

[ref15] 15. Gulbudak H, Martcheva M. Forward hysteresis and backward bifurcation caused by culling in an avian influenza model. Mathematical Biosciences. 2013;246(1):202–212. pmid:24056242
View Article
PubMed/NCBI
Google Scholar

[39] View Article

[40] PubMed/NCBI

[41] Google Scholar

[ref16] 16. Shim E, Galvani AP. Evolutionary repercussions of avian culling on host resistance and influenza virulence. PLoS ONE. 2009;4(5):1–8.
View Article
Google Scholar

[43] View Article

[44] Google Scholar

[ref17] 17. Kim T, Hwang W, Zhang A, Sen S, Ramanathan M. Multi-agent model analysis of the containment strategy for Avian Influenza (AI) in South Korea. Proceedings—IEEE International Conference on Bioinformatics and Biomedicine, BIBM 2008. 2008; p. 335–338.

[ref18] 18. Kim T, Hwang W, Zhang A, Sen S, Ramanathan M. Multi-agent modeling of the South Korean avian influenza epidemic. BMC Infectious Diseases. 2009;10(5):236.
View Article
Google Scholar

[47] View Article

[48] Google Scholar

[ref19] 19. Boender GJ, Hagenaars TJ, Bouma A, Nodelijk G, Elbers ARW, De Jong MCM, et al. Risk maps for the spread of highly pathogenic avian influenza in poultry. PLoS Computational Biology. 2007;3(4):704–712.
View Article
Google Scholar

[50] View Article

[51] Google Scholar

[ref20] 20. Backer JA, van Roermund HJW, Fischer EAJ, van Asseldonk MAPM, Bergevoet RHM. Controlling highly pathogenic avian influenza outbreaks: An epidemiological and economic model analysis. Preventive Veterinary Medicine. 2015;121(1-2):142–150. pmid:26087887
View Article
PubMed/NCBI
Google Scholar

[53] View Article

[54] PubMed/NCBI

[55] Google Scholar

[ref21] 21. Bett B, Henning J, Abdu P, Okike I, Poole J, Young J, et al. Transmission rate and reproductive number of the H5N1 highly pathogenic avian influenza virus during the december 2005-July 2008 epidemic in Nigeria. Transboundary and Emerging Diseases. 2014;61(1):60–68. pmid:22925404
View Article
PubMed/NCBI
Google Scholar

[57] View Article

[58] PubMed/NCBI

[59] Google Scholar

[ref22] 22. Stegeman A, Bouma A, Elbers A, de Jong M, Nodelijk G, de Klerk F, et al. Avian Influenza A Virus (H7N7) Epidemic in The Netherlands in 2003: Course of the Epidemic and Effectiveness of Control Measures. The Journal of Infectious Diseases. 2004;190(12):2088–2095. pmid:15551206
View Article
PubMed/NCBI
Google Scholar

[61] View Article

[62] PubMed/NCBI

[63] Google Scholar

[ref23] 23. Ward MP, Maftei D, Apostu C, Suru A. Estimation of the basic reproductive number (R0) for epidemic, highly pathogenic avian influenza subtype H5N1 spread. Epidemiology and Infection. 2009;137(2):219–226. pmid:18559127
View Article
PubMed/NCBI
Google Scholar

[65] View Article

[66] PubMed/NCBI

[67] Google Scholar

[ref24] 24. Dorigatti I RRPABL Mulatti P. Modelling the spatial spread of H7N1 avian influenza virus among poultry farms in Italy. Epidemics. 2010;2(1):29–35. pmid:21352774
View Article
PubMed/NCBI
Google Scholar

[69] View Article

[70] PubMed/NCBI

[71] Google Scholar

[ref25] 25. Darbon A, Valdano E, Poletto C, Giovannini A, Savini L, Candeloro L, et al. Network-based assessment of the vulnerability of Italian regions to bovine brucellosis. Preventive Veterinary Medicine. 2018;158(March):25–34. pmid:30220393
View Article
PubMed/NCBI
Google Scholar

[73] View Article

[74] PubMed/NCBI

[75] Google Scholar

[ref26] 26. Open Data Porta. Eumsung-gun domestic animals statistical research;. Available from: http://www.data.go.kr.

[ref27] 27. National Geographic Information Institute. Administrative boundaries;. Available from: http://data.nsdi.go.kr/dataset/20180927ds0059.

[ref28] 28. Begon M, Bennett M, Bowers RG, French NP, Hazel SM, Turner J. A clarification of transmission terms in host-microparasite models: numbers, densities and areas. Epidemiology and infection. 2002;129(1):147–153. pmid:12211582
View Article
PubMed/NCBI
Google Scholar

[79] View Article

[80] PubMed/NCBI

[81] Google Scholar

[ref29] 29. McCallum H, Barlow N, Hone J. How should pathogen transission be modeled? Trends in Ecology and Evolution. 2001;16(6):295–300.
View Article
Google Scholar

[83] View Article

[84] Google Scholar

[ref30] 30. Ministry of Agriculture, Food and Rural Affairs. HPAI emergency action guidelines;. Available from: http://www.mafra.go.kr.

[ref31] 31. Heffernan JM, Smith RJ, Wahl LM. Perspectives on the basic reproductive ratio. Journal of the Royal Society, Interface / the Royal Society. 2005;2(4):281–293.
View Article
Google Scholar

[87] View Article

[88] Google Scholar

[ref32] 32. van den Driessche P. Reproduction numbers of infectious disease models. Infectious Disease Modelling. 2017;2(3):288–303. pmid:29928743
View Article
PubMed/NCBI
Google Scholar

[90] View Article

[91] PubMed/NCBI

[92] Google Scholar

[ref33] 33. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences. 2002;180(1):29–48. pmid:12387915
View Article
PubMed/NCBI
Google Scholar

[94] View Article

[95] PubMed/NCBI

[96] Google Scholar

[ref34] 34. Diekmann O, Heesterbeek JAP, Metz J. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology. 1990;28(4):365–382. pmid:2117040
View Article
PubMed/NCBI
Google Scholar

[98] View Article

[99] PubMed/NCBI

[100] Google Scholar

[ref35] 35. Assessing the probability of introduction and spread of avian influenza (AI) virus in commercial Australian poultry operations using an expert opinion elicitation. PLoS ONE. 2018;13(3):1–19.
View Article
Google Scholar

[102] View Article

[103] Google Scholar

[ref36] 36. Gilbert M, Pfeiffer DU. Risk factor modelling of the spatio-temporal patterns of highly pathogenic avian influenza (HPAIV) H5N1: A review. Spatial and Spatio-temporal Epidemiology. 2012;3(3):173–183. pmid:22749203
View Article
PubMed/NCBI
Google Scholar

[105] View Article

[106] PubMed/NCBI

[107] Google Scholar

[ref37] 37. Hill SC, Lee YJ, Song BM, Kang HM, Lee EK, Hanna A, et al. Wild waterfowl migration and domestic duck density shape the epidemiology of highly pathogenic H5N8 influenza in the Republic of Korea. Infection, Genetics and Evolution. 2015;34:267–277. pmid:26079277
View Article
PubMed/NCBI
Google Scholar

[109] View Article

[110] PubMed/NCBI

[111] Google Scholar

[ref38] 38. Manore C, McMahon B, Fair J, Hyman JM, Brown M, Labute M. Disease properties, geography, and mitigation strategies in a simulation spread of rinderpest across the United States. Veterinary Research. 2011;42(1):1–12.
View Article
Google Scholar

[113] View Article

[114] Google Scholar

[ref39] 39. Meadows AJ, Mundt CC, Keeling MJ, Tildesley MJ. Disentangling the influence of livestock vs. farm density on livestock disease epidemics. Ecosphere. 2018;9(7).
View Article
Google Scholar

[116] View Article

[117] Google Scholar

[ref40] 40. Tildesley MJ, Bessell PR, Keeling MJ, Woolhouse MEJ. The role of pre-emptive culling in the control of foot-and-mouth disease. Proceedings Biological sciences / The Royal Society. 2009;276(1671):3239–48.
View Article
Google Scholar

[119] View Article

[120] Google Scholar

[ref41] 41. Backer JA, Hagenaars TJ, Nodelijk G, van Roermund HJW. Vaccination against foot-and-mouth disease I: Epidemiological consequences. Preventive Veterinary Medicine. 2012;107(1-2):27–40. pmid:22749763
View Article
PubMed/NCBI
Google Scholar

[122] View Article

[123] PubMed/NCBI

[124] Google Scholar

[ref42] 42. Buhnerkempe MG, Tildesley MJ, Lindström T, Grear DA, Portacci K, Miller RS, et al. The impact of movements and animal density on continental scale cattle disease outbreaks in the United States. PLoS ONE. 2014;9(3):1–10.
View Article
Google Scholar

[126] View Article

[127] Google Scholar

Figures

Abstract

Introduction

Materials and methods

Epidemiological and geographical data

Culling rate functions

PE culling rate function.

IP culling rate function.

Metapopulation model with two species

Reproductive numbers

Results

Parameter estimations

Reproductive numbers

PE culling strategy

Discussion

Conclusions

Supporting information

S1 Table. Time series outbreak data (chicken farms) of the sixteen towns [1].

S2 Table. Time series outbreak data (duck farms) of the sixteen towns [1].

S3 Table. Pairwise distance matrix of the sixteen towns (km).

S4 Table. Daily AI report and IP culling number.

References