Metapopulation dynamics enable persistence of influenza A, including A/H5N1, in poultry.

Highly pathogenic influenza A/H5N1 has persistently but sporadically caused human illness and death since 1997. Yet it is still unclear how this pathogen is able to persist globally. While wild birds seem to be a genetic reservoir for influenza A, they do not seem to be the main source of human illness. Here, we highlight the role that domestic poultry may play in maintaining A/H5N1 globally, using theoretical models of spatial population structure in poultry populations. We find that a metapopulation of moderately sized poultry flocks can sustain the pathogen in a finite poultry population for over two years. Our results suggest that it is possible that moderately intensive backyard farms could sustain the pathogen indefinitely in real systems. This fits a pattern that has been observed from many empirical systems. Rather than just employing standard culling procedures to control the disease, our model suggests ways that poultry production systems may be modified.

Introduction

Highly pathogenic influenza A/H5N1, colloquially referred to as “Avian Flu”, was first identified in Hong Kong in 1997 as a cause of fatal respiratory illness in humans [1]. After spreading from Asia to Europe and Africa in 2005, it has persisted since then. It has been reported from wild waterfowl in Asia [2] and Europe [3]–[5], as the cause of outbreaks of poultry disease in Asia, Europe and Africa [6]–[12], and as the cause of repeated zoonotic transmission to humans in Asia, Europe and Africa [13]–[16]. Although human-to-human transmission is rare, this strain of avian flu represents a significant pandemic threat, with 360 reported human deaths by the end of 2012 [17], and we now know artificial selection has been able to create a strain of A/H5N1, that has only a modicum of mutations from naturally occurring strains, that can make the virus transmissible among mammals while remaining pathogenic [18]. Yet it remains unclear why this has not occurred outside the laboratory. In addition, the case fatality rates for human A/H5N1 infections remain high [17], as do the mortality rates in chickens, some breeds of domestic ducks, and some species of wild birds [19], making control measures difficult and urgent [20]. Although many agencies and governments have tried to control A/H5N1 spread with some degree of success, it has not been eradicated nor displaced by a less pathogenic strain of Influenza [21], [22].

Indonesia, Egypt and Vietnam represent the three main foci of human A/H5N1, with over 75% of human cases occurring there [17]. Twelve other countries have reported sporadic and occasional human cases and 48 additional countries have reported infected animals [17], [23] with no reported spillover to people. Thus a key question is how A/H5N1 has persisted across the Eastern Hemisphere, but with only very distinct, geographically disparate foci of human cases. One explanation for this may be differences in the nature of poultry farming in affected versus unaffected countries. In particular, a number of studies have proposed backyard poultry rearing as a key risk factor for human infection [8], [24]–[27]. However, many countries that have dense human populations, extensive backyard poultry rearing, and repeated presence of A/H5N1, have reported no or only a few human cases (e.g. Bangladesh, India) [17].

At a large spatial scale, the presence of A/H5N1 in Asia has been correlated with duck density and rice cropping patterns [21], [28], [29], and its spread has been linked to bird migration and global poultry trade [30]–[36]. However, because of high livestock and human mortality, and the short duration of outbreaks, A/H5N1 dynamics have been difficult to study. Here we examine the spatial dynamics of A/H5N1 using theoretical models of disease spread and attempt to understand its ability to persist in different types of domestic poultry rearing operations.

It is a well-established theory in ecology that spatial structuring of a population, as opposed to a single well-mixed population, may alter dynamics [37], [38]. Spatial structure can dampen cycles, allow for co-existence of pathogens with hosts, and generally enable the long-term persistence of often lethal pathogens within a host population [39]–[41]. Employing a metapopulation approach has been successfully used to understand the dynamics of measles in small cities and towns [42], the sporadic nature of Hendra virus outbreaks [43], rabies in wild and domestic dogs [44], and cattle diseases [45]. In the current paper, we apply this understanding specifically to the case of A/H5N1 in domestic ducks, and poultry more broadly. We will demonstrate that spatial structuring can allow A/H5N1 to persist for an indefinite period of time without re-introduction from wild bird populations. We will also investigate the effects of some possible control measures, such as culling and cleaning. While these measures may help reduce the impact and spillover risk of A/H5N1, they may not be able to eliminate unless the effort put into making them effective is enormous.

Methods: Model

We designed our model to represent a single species in several patches (x ∈ {0,…J]}. The model assumes a network of patches, representing farms, markets, or traders, linked together. Our model is a stochastic simulation developed using the Gillespie algorithm [46], [47] because we are focused on variability in persistence, but we present the deterministic skeleton of the model for simplicity and clarity. Within a patch, the model uses standard SIR dynamics for each species, with the addition of the presence of a common environmental reservoir of infection [48]. Within each patch (x), the dynamics follow these equations:

Where S represents the number of susceptibles, I infected, R recovered, and V virions in the environment. The initial conditions for all simulations are an entirely susceptible population, except for two infected individuals, located in two randomly chosen patches. We have included the indirect transmission model of Breban et al. [48], where each infected host sheds virus into the environment at rate (), and virions in the environment degrade at rate (η). We use only the only the warm temperature parameter values (see Table 1), as we are focused here on poultry in subtropical and tropical environments. Transmission may occur indirectly from the environment , with a dose-dependence of likelihood of infection , which depends on virion infectiousness (). Direct transmission occurs at rate β. Infected animals recover at rate (), or die due to the disease at rate (α). In addition, there may be non-influenza related mortality (μ), for all categories (Sx,Ix,Rx) of animals. We do not include birth, because for most real poultry systems, chicks are introduced by other means. Although not shown, the model keeps track of the number of dead animals from each infection category. For the introduction of a single infected individual into a fully susceptible population (S0), without any infectious virions in the environment, this model has a within-patch basic reproduction number of:

Table 1

Parameters of the model.

Symbol	Value/Range	Definition	Units	Reference
	0.004	direct transmisison	animal−1 day−1	[73]
	0–0.167	recovery rate	day−1	[51]–[53], [73]
	0–0.4	disease mortality rate	day−1	[51]–[53], [73]
	0.001	environmental uptake rate	day−1	[48]
	1.96 ⋅ 10−4	virion infectiousness	virion−1	Calculation from [48], ID50 from [74]
	0.14	virus degradation rate	day−1	[48]
	105	shedding rate	day−1	[48]
	0.00–0.30	movement rate	day−1	[59]
	2	network neighborhood	patches	–
	0.0596	re-wiring probability		–

Animals may leave the patch at rate (ω) from any class, presumably by human activity in the case of poultry. Animals that leave their current patch are then moved to another patch according to the function:

Where χ(x,y) represents the probability of movement between patches x and y, given a departure from patch y, and N represents the appropriate population (susceptible, infected or recovered), with the constraint that all χ(x,y) must sum to 1 across all possible source patches. Here we use the row normalized adjacency matrix of a network model to specify the interconnectedness of the system. In the absence of relevant empirical data, we have examined two classic random networks. We have used the Watts-Strogatz random network model, known as the small-world model with neighborhood and re-wiring probability for global connections (e.g., Figure 1a) [49]. In order to examine the impact of network structure, we have also used a Barabasi random network (e.g., Figure 1b) [50].

Figure 1

Two example random networks.

(A) A Watt-Strogatz random small world network with 64 vertices, with  = 2.33 and ρ = 0.0596 (as in Table 1), (B) A Barabasi random network with 64 vertices, with the power of the preferential attachment set to 0.5 and the zero appeal to 1.0, both represent a network of 64 farms and markets, each with 500 animals.

Additionally, we have implemented a model of passive surveillance based control in the model. If the number of dead birds within a given period (τCrit) exceeds a threshold (ICrit), it is presumed an authority would be notified given a probability of reporting (πReport), and that authority would test the living animals for infection, with a given probability of detection (πDetect). If the infection is detected, then all animals can be culled, and/or the environment cleaned ( set to 0).

Although we have developed a full multispecies model, in order to focus on the effects of population structure, we limit our consideration to a single species of farmed livestock, namely ducks [51], [52]. We do examine different parameter values, which cover the case of chickens as well [53]. We also treat the total population size as fixed at first. Population size becomes variable as animals move through the network, increasing as animals move in from other patches, as well as decreasing due to mortality and animals moving out of patches. We did this so as not to confound the results of spatial structure with those of population size. We generally examine equal initial population sizes, again to focus our results on structure; but we have looked a few key cases with mixed population sizes.

Results

Our key finding (Figure 2) is that persistence of A/H5N1 is highly dependent on farm size. For very small patches, the within-patch R0 of the pathogen is substantially less than one and most epidemics fail, i.e. there are no cases of secondary transmission. This is true across different population structures (Figures S1–S5). The qualitative result is also robust to different parameter values (Figures S6–S9) [51]–[53]. Although truly random networks can require very large patches for persistence, this may be because the initial two infected hosts frequently land in patches that are completely disconnected from the rest of the network, and it is only once the network is small enough that the infection can actually spread throughout the network (Figure S2). In contrast, the transitions are very stark for a pure nearest neighbor network, the other extreme of the Watts-Strogatz model (Figure S1). The results are qualitatively robust for the topologically very different Barabasi network (Figure S3), but this does not alter the general qualitative conclusions. Including a few larger farms in a matrix of small farms does not alter the results very much until there is a gigantic farm, which reduces the variability in outcomes (Figure S4 & S5).

Figure 2

For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations.

Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program. Gray area represents parameter region where 1

At low patch sizes, if there are some cases of secondary transmission, the epidemic dies out rapidly due to stochasticity (e.g., Figure 3a). For very large patches, the within-patch R0 of the pathogen is so high that the epidemic rapidly burns through the population, and as long as the network is sufficiently connected that all patches become infected, they all experience a full epidemic (e.g., Figure 3c). In instances where patches are of moderate size however, resulting in R0 values between 1 and approximately 6, patches experience asynchronous local mini-epidemics (e.g., Figure 3b), which prolong the global epidemic. However in these cases, the epidemic doesn’t necessarily infect a large portion of the population, and if pathogenicity is low, may not cause noticeable outbreaks. Variation in the rates of direct transmission, mortality, and recovery, as well as the details of environmental transmission, can alter the pathogen’s ideal patch size, but does not change the fundamental result that a metapopulation of moderate sized farms is best for persistence. Simple mixed initial population size models, with either one (Figure S4) or five (Figure S5) larger populations, and several smaller moderately sized populations (N = 250) do not alter the conclusions until these larger populations are quite large (>5000) relative to the total population size. Even then, the effect is more to reduce the variability, rather than a strong alteration of the mean. This highlights the role that the moderate size farms have in allowing the epidemic to smolder, even if there are a few large epidemics in the system without any control measures.

Figure 3

Dynamics over time of for three different simulation runs, solid grey line represents global prevalence of infection across all patches (farms and markets), colored dashed lines represent abundance of infected individuals within each patch (farm and market).

(A) 1280 patches of size 25, longest simulation run, note only three patches infected, (B) 128 patches of size 250, random simulation run, note global prevalence always less than 2.5%, (C) 16 patches of size 2000, random simulation run, note near deterministic similarity of epidemic in each patch.

In Figure 4, we implement a control program, namely a passive surveillance program that culls animals and cleans the environment when influenza is detected. While these measures can substantially reduce the total number of individuals infected in larger patches, they are ineffective at curtailing persistence in a network of moderate sized populations. Figure S9 demonstrates that our conclusions are general, although the quantitative details depend on exact parameters.

Figure 4

For a fixed total single species total population size (32,000), without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations.

Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Here with control measures implemented, (πReport = 0.9, πDetect = 0.9, τCrit = 1, ICrit = 5). Gray area represents parameter region where 1

Discussion

Our model demonstrates a strong role of farm size in the risk of highly pathogenic avian influenza (HPAI) outbreaks. Our results suggest that true small-scale subsistence farming at low densities has a low risk of HPAI outbreaks, and since the risk of outbreaks is reduced, the risk of spillover to people should be reduced. These results run counter to a number of previous studies which have linked human cases of highly pathogenic A/H5N1 influenza to small-scale, backyard poultry operations and proposed that this is a major risk factor for future spillover events [54]–[57], however, the small and moderate scales (discussed below) in our model are rarely separated in the empirical literature, as both fall into the FAO/OIE sector 4 designation [58], [59]. In most cases, these studies consist of investigations of human cases, outbreaks, or policy recommendations based on epidemiological analyses of these. Thus results here suggest that patchily distributed very small-scale low density poultry production are insufficient to sustain epidemics, and a fragmented trade network may instead reduce the probability of sustained transmission.

As human density increases, patches of poultry likely become larger – either from the development of modest poultry production by local entrepreneurs, or because individual family flocks intermingle so extensively to become a single flock in high-density environments, or from a mix of these scenarios. In this case, a smoldering epidemic might allow an HPAI to persist without ever causing enough livestock mortality to enable effective intervention. The countries that have had the most reported H5N1 cases share these characteristics in common, despite geographic distance. Indonesia, Vietnam, and Egypt have high human population density, where people raise poultry both for subsistence and for income, but on a moderately small, household scale [58].

Our modeling work suggests that truly large-scale poultry production, such as Concentrated Animal Feeding Operations (CAFOs), could reduce the persistence of A/H5N1 by reducing the likelihood of undetected and under-reported smoldering epidemics and thus reduce the risk of spillover to humans. A critical part of this would be a control program that succeeds at identifying outbreaks in these farms sufficiently early to provide a public health benefit; otherwise smoldering persistence would become raging epidemics. Additionally, it would need to be clear that the advantages from making epidemics more detectable offset the disadvantage of any increase in the size of epidemics for the risk of spillover to humans. A small number of studies have highlighted the lack of biosecurity in intensive poultry farms [60], [61]. However, these studies are focused on developed countries, where there is an expectation of reliable, advanced biosecurity, and when lapses are identified, they are considered significant. In developing countries, the difference in biosecurity between industrial farms and backyard production is likely far larger, even if neither is as biosecure as in developed countries. Here, there are potential practical reasons for a protective role of large-scale farms. In these developing country large-scale farms, there is likely to be less trading in and out of the population and far more surveillance for mortality than in small backyard flocks. Therefore an HPAI outbreak would be large and noticeable and therefore more likely to be subject to control measures such as the de-flocking or mass vaccination strategies used widely in China, southeast Asia and Egypt [58].

Comparing the results from different network structures, it is clear that this is an important factor, particularly for systems of larger farms. A nearest neighbor network (no long-distance links across the network) produced results almost similar to culling. However, increasing farm size did reduce persistence in all cases, albeit with different patterns. Yet the critical threshold farm-size for enabling persistence was similar across network patterns, demonstrating that the within patch R0 is the most important factor for this threshold. Unfortunately, the detailed comparative data on real trade networks, and poultry farming practices across countries that would be needed to test this model are currently lacking, or relatively inaccessible proprietary data, although we are seeking to collaborate with other groups and nations to obtain this data, as well as working on our own empirical research.

It is important to note that our model can explain the persistence of A/H5N1 influenza without invoking repeated introductions from wildlife, or domestic, reservoirs. Wild waterfowl harbor a diversity of influenza strains, including A/H5N1, which has been reported from wild birds in Egypt [62], South Asia [32], China [2], [34], [63], Europe [64] and Africa [65], [66]. However, while there is evidence for a role of wild birds in the spread and introduction of A/H5N1 [31], it is not clear what role infected wild birds play in persistence of HPAIs generally. Our modeling suggests that the wild birds are not required for persistence, in line with a recent empirical review by Gauthier-Clerc [67].

Culling remains the most widespread and commonly used approach for dealing with HPAI infections in endemic countries [58]. Our model demonstrates the effectiveness of culling in reducing the number of infected individuals in large poultry populations. However, neither culling, nor culling and cleaning of the environment was able to reduce persistence of influenza in our simulations for metapopulations of moderate sized farms. Our model results suggest that changes to the type of farms present within endemic countries would have a more significant impact on the persistence of HPAI, and therefore the long-term effectiveness of control programs. In essence, once the protein demands of country require the intensification of poultry production beyond subsistence, from an emerging disease risk prospective it may be best to move to CAFOs as quickly as possible. A preliminary suggestion in terms of network structure would be to subdivide the farm to market chain into as many small separate networks as possible, but this may not be feasible.

In this study, we did not examine the efficacy of vaccination as a control strategy because of the complexity of escape mutations, and problems with its long-term use as a control measure [20], [68], [69]. Therefore, a model that proposed to examine a vaccination control strategy would have to include strain variation and antigenic selection (e.g. [70]), which is beyond the scope of this work. Future work should investigate the ability of this model to fit observed dynamics in empirical systems, the potential effects of variability of initial farm size, and alternate network structure [71], both empirical and theoretical. This model could be helpful in optimizing control efforts could potentially yield tremendous benefits to the poultry industry, as the virus depresses exports in affected countries, with China and Thailand alone losing $900 million due A/H5N1 from to 2003 to 2005 [72].

Our study demonstrates that moderate size patches of poultry may substantially contribute to the persistence of influenza in countries where such production is a dominant force in the poultry production system. Although passive surveillance and culling may reduce prevalence, and infection burden, and thus risk of spillover it is unlikely to eliminate influenza from these countries. Development of more bio-secure, intensively monitored, larger scale commercial poultry production may be the best route to risk reduction in countries with substantial need for intensive poultry production.

Alternate Watts-Strogatz network, here re-wiring probability ρ = 0, thus only the two nearest neighbors in either direction are connected, and there are no random long distance connections across the network. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program. Gray area represents parameter region where 1< R0<6 within a patch.

(TIF)

Click here for additional data file.

Alternate Watts-Strogatz network, here re-wiring probability ρ = 0, which transforms the network to an essentially randomly wired network. This can create sub-networks that are disconnected from the majority of the network leading to high variability in simulation results, and pushing the persistence area to very large farm sizes. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program. Gray area represents parameter region where 1< R0<6 within a patch.

(TIF)

Click here for additional data file.

Alternate random network, here a non-directional Barabasi network [ , with a zero appeal of 1 and a power of preferential attachment of 0.5. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program. Gray area represents parameter region where 1< R0<6 within a patch.

(TIF)

Click here for additional data file.

Standard small world network (ρ = 0.6), but there is now one larger farm of variable size. The top row of the x-axis is the size of the larger patch, while the bottom row is the number of smaller patches, all of size 250 hosts. Thus the far left of the graph is the same as for Figure 1 in the main text. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program.

(TIF)

Click here for additional data file.

Standard small world network (ρ = 0.6), but there is now five larger farms of variable size. The top row of the x-axis is the size of each of the five larger patches, while the bottom row is the number of smaller patches, all of size 250 hosts. Thus the far left of the graph is the same as for Figure 1 in the main text. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except α = 0.1111, and ω = 0.03, without any infection control program.

(TIF)

Click here for additional data file.

Alternate transmission parameters that more closely resemble H5N1 infections in chickens with no recovery and faster mortality [ . For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except β = 0.0081, γ = 0, α = 0.32, and ω = 0.05, without any infection control program.

(TIF)

Click here for additional data file.

Alternate transmission parameters that emphasize environmental transmission more and direct transmission less, with higher mortality. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, including environmental transmission, except β = 0.003, υ = 0.002, α = 0.222, and ω = 0.143, without any infection control program.

(TIF)

Click here for additional data file.

Alternate transmission parameters without environmental transmission. For a fixed total single species total population size 32,000, without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Other parameters are as in Table 1, excluding environmental transmission, i.e. υ = 0.0, without any infection control program.

(TIF)

Click here for additional data file.

Less effective control program, π For a fixed total single species total population size (32,000), without non-influenza mortality (μ = 0), the effect of changing local patch size and patch number on (A) frequency of epidemic failure, (B) median length of epidemic in days, and (C) median total number of animals infected over 100 simulations. Dotted lines represent the empirical 97.5% and 2.5% percentiles, creating a 95% bootstrap confidence interval. Here with control measures implemented, (πReport = 0.1, πDetect = 0.9, τCrit = 1, ICrit = 5).

(TIF)

Click here for additional data file.