Numerous empirical and modeling studies dealt with the quantification of indirect, particularly airborne transmission of FMD and other diseases between domestic livestock holdings [e.g., Ref. (71, 80–84), reviewed for FMD in Ref. (85)], but IT of FMDV in wildlife animals has, to our knowledge, not yet been quantified. We developed a modeling approach that breaks down IT into components that are accessible to experimental measurements, namely pathogen shedding, survival/decay in the environment, contact with infectious material, and infection according to a dose–response relationship. Although some experiments quantified pathogen excretion and secretion of FMDV [reviewed in Ref. (12)] and other pathogens [see, e.g., Ref. (86) for CSF] by domestic animals, knowledge for wildlife is rare (87). The large differences between domestic animal species regarding the shedding rates of FMDV (12, 75) call for further attention to this issue. The same applies for the susceptibility of different species, i.e., the dose–response relation (12, 75). Some quantification for domestic animals can be found in the literature [e.g., Ref. (70, 88) for FMD], but the qualitative relation between dose in the environment and probability of infection is often unclear (89, 90). Survival outside the host has been investigated for several pathogens in animal products and excrements under laboratory conditions (for references, see INTRODUCTION), but further research is necessary for environmental factors that influence pathogen survival. The contact of animals with viral contamination in the environment remains the most uncertain parameter. Here, an inverse parameter fitting approach could aid the quantification. Given assumptions for the other parameters, contact to viral contamination could be estimated from the probability of infection.

Experimental investigations of virus survival outside the host depict striking dependence on temperature and humidity [see, e.g., Ref. (26) for FMD, CSF, BVDV, and PPV; (91) for PRRS virus; (92) for influenza A]. This fact gives rise to seasonal fluctuations of

#### **REFERENCES**

the magnitude of IT. Indeed, for several viral diseases, fluctuations of their transmission were associated with climatic seasonality, partly related to virus survival outside the host [see, e.g., Ref. (92, 93) for influenza viruses; (94) for hepatitis A]. Therefore, climatic factors are expected to play a role in regional variations of the epidemiology of infectious diseases with an IT mode.

With this work, we contribute to the research on IT, which is still in an early stage but attracting increasing attention. Previous work focused on the impact of IT on key figures of host–pathogen systems such as the basic reproductive number (20), disease persistence (41), and formal conditions of relevance for modeling (40).

Our results resemble findings byWearing et al. (1) and Almberg et al. (20), which show that a neglect of prolonged infectiousness, e.g., through environmental pathogen stages or inappropriate assumptions about the infectious period, may result in an underestimate of *R*0, if derived from the prevalence growth rate. Reciprocal, in our study prevalence growth rates decreased under IT despite equal reproductive potential (*R*0). Thus, we transferred the findings regarding the relevance of IT from a theoretical underestimation of infection dynamics, i.e., *R*0, to the application-oriented context of designing surveillance of any particular wildlife disease, i.e., *R*0 being fixed.

We conclude that a simplified aggregation of transmission processes, particularly a neglect of environmental pathogen stages, may considerably bias model predictions of the performance of disease surveillance and mitigation strategies. We state that this applies even for pathogens with an average environmental survival time that is comparatively short compared to the infectious period of the host.

#### **AUTHOR CONTRIBUTIONS**

ML and H-HT conceived and designed the experiments; ML performed the experiments and analyzed the data; and ML, SK-S, and H-HT developed the model and wrote the manuscript.

#### **FUNDING**

ML was partially funded by the European Food Safety Agency (EFSA).

#### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fvets.2016. 00110/full#supplementary-material.

Islands: implications for transmission to wildlife. *Auk* (2008) 125(2):445–55. doi:10.1525/auk.2008.06235

swine fever – legend or actual epidemiological process? *Prev Vet Med* (2012) 106:185–95. doi:10.1016/j.prevetmed.2012.01.024

porcine reproductive and respiratory syndrome virus in aerosols. *Vet Res* (2007) 38:81–93. doi:10.1051/vetres:2006044

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2016 Lange, Kramer-Schadt and Thulke. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Impact of Network Activity on the Spread of Infectious Diseases through the German Pig Trade Network

#### *Karin Lebl1 \*, Hartmut H. K. Lentz1 , Beate Pinior2 and Thomas Selhorst3*

*1 Institute of Epidemiology, Friedrich-Loeffler-Institute, Greifswald, Insel Riems, Germany, 2 Institute for Veterinary Public Health, University of Veterinary Medicine Vienna, Vienna, Austria, 3Unit Epidemiology, Statistics and Mathematical Modelling, Federal Institute for Risk Assessment, Berlin, Germany*

#### *Edited by:*

*Tariq Halasa, Technical University of Denmark, Denmark*

#### *Reviewed by:*

*Lina Mur, Kansas State University, USA Beatriz Martínez-López, UC Davis, USA*

> *\*Correspondence: Karin Lebl<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Veterinary Epidemiology and Economics, a section of the journal Frontiers in Veterinary Science*

*Received: 17 March 2016 Accepted: 07 June 2016 Published: 21 June 2016*

#### *Citation:*

*Lebl K, Lentz HHK, Pinior B and Selhorst T (2016) Impact of Network Activity on the Spread of Infectious Diseases through the German Pig Trade Network. Front. Vet. Sci. 3:48. doi: 10.3389/fvets.2016.00048*

The trade of livestock is an important and growing economic sector, but it is also a major factor in the spread of diseases. The spreading of diseases in a trade network is likely to be influenced by how often existing trade connections are active. The activity α is defined as the mean frequency of occurrences of existing trade links, thus 0 < α ≤ 1. The observed German pig trade network had an activity of α = 0.11, thus each existing trade connection between two farms was, on average, active at about 10% of the time during the observation period 2008–2009. The aim of this study is to analyze how changes in the *activity* level of the German pig trade network influence the probability of disease outbreaks, size, and duration of epidemics for different disease transmission probabilities. Thus, we want to investigate the question, whether it makes a difference for a hypothetical spread of an animal disease to transport many animals at the same time or few animals at many times. A SIR model was used to simulate the spread of a disease within the German pig trade network. Our results show that for transmission probabilities <1, the outbreak probability increases in the case of a decreased frequency of animal transports, peaking range of α from 0.05 to 0.1. However, for the final outbreak size, we find that a threshold exists such that finite outbreaks occur only above a critical value of α, which is ~0.1, and therefore in proximity of the observed activity level. Thus, although the outbreak probability increased when decreasing α, these outbreaks affect only a small number of farms. The duration of the epidemic peaks at an activity level in the range of α = 0.2–0.3. Additionally, the results of our simulations show that even small changes in the activity level of the German pig trade network would have dramatic effects on outbreak probability, outbreak size, and epidemic duration. Thus, we can conclude and recommend that the network activity is an important aspect, which should be taken into account when modeling the spread of diseases within trade networks.

Keywords: network analysis, disease spread, trade activities, temporal network, animal movements, epidemiology

#### INTRODUCTION

Live animal trade represents an important economic sector but is permanently subject to fluctuations. For instance, consignments of pigs increased to 48% within EU-27 member states between 2005 and 2009 (1). However, the financial crisis in the subsequent years might have lessened this effect. The importance of live animal trade on the economy is also demonstrated during animal disease outbreaks. Trade restrictions with movement bans cause enormous financial losses for the affected livestock holdings and countries. For example, the outbreak of classical swine fever (CSF) in the 1990s in Germany led to an economical loss of approximately €1 billion (2). Thus, as demonstrated during CSF outbreak in Germany, livestock trade between farms is one of the major routes for the spread of animal diseases, although other infection routes, like proximity to infected herds or contact with contaminated persons and vehicles, exist as well (2).

Scientific research has primarily focused on the influence of the trade structure of farms on disease dynamics (3, 4). Farms differ with respect to their trade activity, i.e., with respect to the number of trading partners, trade connections, trade volume, and time intervals (5). Within the trade network, farms with greater trade activities are the most important contributors to disease spread (6). Veterinary epidemiology assessments utilized social network analysis (SNA) tools, such as centrality measures, developed within the field of social sciences, to calculate the importance of farms for the spread of animal infections. Numerous centrality measures, such as in- and out-degree, betweenness, and closeness (7), were correlated with standard epidemiological parameters, such as size of an epidemic, duration of the epidemic, time to peak of the epidemic, and the basic reproduction number *R*0 (4, 8–10).

Previous studies applying SNA on pig trade networks have already provided important insight for disease prevention and control. One aspect of this research was the identification of the structure of trade communities (11, 12). Another essential finding was that there is a large degree of heterogeneity associated with movements of pigs at the movement level and at the premise specific network level as well (13). As a result, pig trade has a right-skewed distribution of all centrality parameters, i.e., few holdings have high centrality, while most have a low centrality. Thus, strategic removal of the most central nodes would result in a decomposition of the network into fragments, which would interrupt infection chains and prevent further disease spread (14–16). It was also shown that the holding types differ in their centrality measures, which allow for a targeted removal of specific holding types in the case of a disease outbreak (16–18). Further, SNA has been utilized to simulate the spread of specific diseases to estimate the effects of an outbreak, e.g., the spread of Methicillin-resistant *Staphylococcus aureus* (MRSA) through the Danish pig trade network (19).

Although SNA provides useful insights into epidemic dynamics on trade systems, the methods used in SNA do not take into account the temporal ordering of trade links. Whenever a network is traversed using trade links, each traversal has to follow a causal sequence of connections. This constraint can have a significant impact on the spreading paths for pathogens in networks (20). For this reason, recent work has been focused on *temporal network analysis*, where each connection has a time stamp marking its occurrence time. The probability of contagion between two individuals is not constant in time and depends, beside the transmission rate and infectious period, also on the frequency and duration of the contact (21–24). Studies that considered the heterogeneity and duration of contacts and their importance for the epidemic showed the importance to elucidate the time dependency of activities in order to investigate disease dynamics (22, 24–26). Previously, it has been shown that the aggregation of trade links into static networks leads to an overestimation of the epidemic size (27–30), the outbreak probability (31), and the epidemic duration. Thus, scientific research in the veterinary field has increasingly focused on time-dependent networks. Methods have been adapted and extended from static analyses to time-dependent analyses (20, 27, 31–37).

A temporal network view on livestock trade networks includes the frequency of trade links. For the whole system, this frequency can be considered as the pace of trading. This raises the question, whether it makes a difference for a potential spread of an animal disease to transport many animals at the same time (low frequency) or few animals at many times (high frequency). From the economic point of view, it is appropriate to choose a low trade frequency and transport many animals at the same time.

In this work, we analyze the impact of the overall trade frequency on the spread of infectious disease. Hereby, we keep the total trade volume of the network constant and systematically investigate the impact of a changing frequency of traded animals. We define the *activity* of a network by averaging the frequency of all existing trade connections between node pairs and analyze how changes in the activity influence the probability of a disease outbreak, the final outbreak size, and the duration of an epidemic. A discrete stochastic SIR model is used to simulate the spread of a hypothetical disease through the trade network of the German pig production chain.

#### MATERIALS AND METHODS

In order to analyze the influence of network activity on the course of an epidemic, an outbreak model predicting the course of a hypothetical animal disease on a contact network between holdings belonging to the German pig production chain was set up. Besides the outbreak model, we propose a method how to systematically adjust the activity of the network.

#### Data and Network Setup

According to the EU directive EC/2000/15 (38), EU member states are obliged to collect and record livestock movement data in a national database. Pursuant to the German Animal Movement Directive (Viehverkehrsverordnung), each holding in the pig production chain (including piglet production, breeding, raising, fattening, slaughtering, and trading) is obliged to notify the movement of pigs within 7 days. All data are stored in a database, "Herkunftssicherungs- und Informationssystem für Tiere" (HI-Tier). In Germany, movement data for pigs are collected on a daily basis. In general, movement data of livestock comprise information about the source and target farms (unique identifiers), the date of movement, and the number of animals moved (batch size).

For this study, pig movement data from the federal states of Bavaria and Baden-Württemberg between the years 2008 and 2009 were used. It has previously been shown that a period of 2 years is suitable to cover all characteristic properties of the German pig trade network (31). In our data set in most cases (90%), only one movement per week took place between a supplier and buyer. Consequently, we decided to use a weekly timescale for our analysis. In the case of two movements per week, those were merged into one occasion.

To describe the pattern of trade activity over time, a temporal network was constructed. By implementing a temporal network, it is possible to take into account causality for network transversal. In other words, consecutive trade connections have to be temporally ordered in order to make up a valid indirect connection between farms (**Figure 1**). The network comprised nodes and edges, where each edge connected a node pair. Farms were represented by nodes, and movements of animals between farms at a certain point in time were represented by directed edges. A temporal network is defined as (*V*, , *T*), where *V* is the set of nodes within the network, is the set of directed edges, and *T* represents the length of the observation period, as we considered weekly time steps, *T* = 104 weeks. An edge (*u*, *v*, *t*, *w*)∈ describes the movement of *w* pigs from farm *u* to farm *v* at time *t* ≤ *T*. This network comprises |V| = 45,065 and || = 1,237,753 edges (i.e., overall number of transports during the observation period). Further, the static representation of the network was constructed by summing all observations in the temporal network over the study period, such as the static network is the time-aggregated network of . In the static representation of the network *G*(*V*, *E*), *V* represents the set of nodes and *E* the set of directed edges (|*E*| = 112,826). A directed edge between two nodes exists in the static network if a certain animal movement has taken place at least once during the observation period.

The aim of this analysis was to investigate the influence of the network activity on the outbreak size of an epidemic. However, this outbreak size would be strongly influenced by differences in the reachability of the nodes, i.e., nodes form distinct reachability

classes where a significant number of nodes may only cause trivial outbreak (11, 12). To reduce this bias, the data were first tailored to include only nodes, which are, in the static representation of the network, reachable from each other. We used the static network to identify the *largest strongly connected component* (LSCC; in a strongly connected component, each node is reachable by any other node in the component). The further analysis was limited to this LSCC, which we denote as *G*\*. Thus, the static representation of the network enables the disease to reach all nodes in finite time, no matter which node is the source of infection. All nodes and edges, which were not elements of the LSCC in the static network, were removed, as well as the corresponding elements from the temporal network. We hereby implicitly assumed that the concept of connectivity (35, 36) is preserved for the temporal network. In the resulting network, pigs moved between |*V*| = 7,455 farms (number of nodes in the LSCC) and |*E*| = 27,149 transport routes (number of edges in the LSCC) were recorded during the observational period, corresponding to || = 315,481 transports in the temporal network.

#### Setting the Network Activity **α**

Starting from the network generated as described above, we changed the activity systematically. The *activity of a single edge* in a temporal network can be described by its frequency, i.e., how often a certain edge was active during the study period divided by the length of the study period. The *network activity* α was defined as the mean edge frequency of a network, with 0 < α ≤ 1. The network activity α of a temporal network (*V*, , *T*) and its according static representation *G*(*V*, *E*) can be calculated as follows.

$$\alpha = \frac{|\mathcal{E}|}{|E| \times T},\tag{1}$$

where || is the number of edges in the temporal network, |*E*| is the number of edges in the aggregated network, and *T* is the observation period.

In order to investigate the influence of α on disease dynamics, we propose a method to systematically change the network activity. Since the results for a network with shifted α should be comparable to the original network, following constraints had to be considered: (i) the aggregated network *G* remained the same for all α, (ii) the total trade volume remained constant for all α, (iii) the temporal sequence of existing trade routes had to be preserved (see details below), and (iv) the observation period *T* was preserved.

In order to highlight the activity of a temporal network, we computed the activity according to Eq. 1 and denoted a temporal network with a certain network activity as α. For our observed network, we found α = 0.11, and we denote the observed network as α=0.11≡\*.

In order to create networks with a reduced α, randomly chosen edges from \*(*V*, , *T*) were removed. According to constraint (i), edges were removed in a way such that each edge of the aggregated network appeared at least once in the newly generated temporal network.

In order to increase α, we first considered our temporal network as a sequence of static network snapshots. In other

*t* = 1, thus before the disease has reached node *x* at *t* = 2.

words, a temporal network consists of an ordered sequence α(*V*, , *T*) = *G*1, *G*2, … , *G*T, where each *G*t∈ is a static snapshot of the temporal network at time *t*. In order to increase α, each snapshot was first duplicated (once or multiple times) and time-shifted by a certain value chosen at random. Second, these snapshots were merged into a new temporal network. In the case of overlapping edges occurring between the same node pairs (i.e., multiple occurrences of directed edges active at the same time; regardless of their edge weights), the edge weights *w* (i.e., number of transported pigs) were averaged. Using this approach, the existing trade routes remain preserved as required by constraint (iii).

In order to satisfy constraint (iv), we used periodic boundary conditions, i.e., for each edge *(u*, *v*, *t*, *w)* = (*u*, *v*, *t* + *T*, *w*). In other words, if the new times exceeded the observation period *T*, the times were shifted by subtracting *T*.

The procedures described above would already be sufficient to change the activity α of the observed network \*. Nevertheless, both procedures would violate constraint (ii), as the overall sum of edge weights changes as well. Therefore, the new edge weights had to be adjusted. During the observation period, a total of *W* = 24,995,162 transported pigs were recorded. The new edge weights for α were normalized, so that the sum of the new edge weights equaled the total of the observed edge weights *W*. Finally, edge weights for the generated network were rounded to a whole number, with the minimum number of pigs per transport set to one [constraint (i)].

For example, in a first step, we duplicated the graph \* once and conducted a 52-week shift (i.e., a shift of 1 year in the duplicate). Thus, an edge active in the original graph at weeks 2, 40, 63, and 92 would be active in a 52-week-shifted graph at weeks 54, 92, 11, and 40. Merging the original with the time-shifted graph would thus result in a graph where this certain edge is active at weeks 2, 11, 40, 54, 63, and 92 (see **Figure 2** for a more detailed example).

Overall, 22 different networks were generated with different activity values, including the original network \* = 0.11. The considered values for α were approximately evenly distributed in the interval (0; 1].

#### Disease Dynamics SIR Model

In order to analyze the influence of network activity α on the course of an epidemic, we simulated the spread of a disease on different temporal networks α with parameter α. Disease dynamics were modeled by applying a stochastic discrete-time SIR model (29). Farms were treated as epidemiological units that are assigned to one of the three epidemiological states: susceptible (*S*), infected (*I*), and recovered (*R*). The infection spread along an

FIGURE 2 | Example for creating a graph with an increased **α** using a time shift of two time steps. In this example, the original graph has |*V*| = 4 nodes (A, B, C, D) and || = 6 directed edges (with *u* as the starting node and *v* as the receiving node), corresponding to |*E*| = 5 in the time-aggregated network. The edges are active at times *t* ϵ {1, 2, … , 5} (numbers next to the drawn edges), thus *T* = 5. The line widths of the edges correspond to the edge weights *w*. Overlapping edges (i.e., edges with the same *u*, *v*, *t*) are marked in red. The newly generated graph has the same number of nodes, but an increased number of edges (|| = 11). Note that, due to rounding errors, the sum of edge weights in the original and the new graph are only approximately equal.

edge (*u*, *v*, *t*, *w*), if at time *t* the supplying node *u* was in state *I*, and the receiving node *v* in state *S*. Thus, a receiving farm could only became infected, if a transport took place from the supplying farm to the receiving farm during the time period in which the supplying farm was in the *I* state. Infectious nodes stayed in the *I* state for μ time steps, thereafter they passed to the *R* state. Nodes in the *R* state remained in this state until the end of a simulation run. Infectious farms infect susceptible farms with probability *pe*.

Due to the fact that certain information was not available, the following model assumptions were made. (I) farms representing the nodes within the network were all treated identically (8, 29, 31, 39). Thus, in this model, the number of animals on the farm, breed, farm type, or farm practices did not have an effect on the transmission dynamics. (II) the epidemiological status does not alter the trade contact structure. The latter is a strong assumption, but it allowed an examination of the influence of network topology on unmanaged disease dynamics (29).

#### Model Parameters

In order to compute the transmission probability *pe* for each edge, we first considered the risk of infection for each transported animal. For every transport from an infected to a susceptible farm, each transported animal has a probability *p* to infect the receiving node. In this work, the probabilities *p* = {0.25, 0.5, 0.75, 1} were considered. The receiving node became infected, when at least one transported animal spread the disease. The probability *pe* can be described with a binomial function *B*(*w*, *p*), whereas the function depends on the parameters edge weight *w* and an animals' transmission probability *p*.

$$\mathcal{P}\_c = P\left(X > 0\right) \sim B\left(\mathcal{W}, \mathcal{P}\right) = 1 - \left(1 - \mathcal{P}\right)^{\text{w}}.\tag{2}$$

A transmission probability of *p* = 1 corresponds to a highly infectious disease: the supplied farm always became infected, independent of the batch size *w*. This corresponds to a worst-case scenario and is therefore often used in studies investigating the spread of diseases within the trade networks (6, 8, 31).

Nodes remain in the *I* state for the *infectious period* μ and then pass to the *R* state. Nodes in the *R* state remained in this state. In this paper, we considered a constant infectious period of μ = 4 weeks [as estimated for cases, such as CSF, African swine fever, foot-and-mouth disease; Ref. (40)].

#### Initial Conditions

In the analysis presented here, the model predicted the disease dynamics for discrete intervals of 1 week. Initially, all farms were in the susceptible state (*S*). At a randomly chosen time, the state of one randomly selected farm was set to infected (*I*).

The disease dynamics were simulated on a temporal network α. All possible start times and initially infected nodes (index nodes) had the same selection probability. The start times were selected from the interval [1; *T* − 40] to avoid that the durations of the epidemics exceed the observation period of 104 weeks. We chose 40 weeks arbitrarily, as the first test runs showed that the duration of the epidemic only rarely exceeded this time period. However, in some cases, the duration of the epidemic still exceeded the study period – those cases were excluded from the further analysis. The simulation stopped when the number of infectious nodes reached 0.

#### Summary of Parameters

For each value of the activity parameter coming from one of the 22 investigated α, each with the 4 transmission probabilities as described above, the simulation was repeated 2,000 times, as test runs showed that this number of iterations provided robust results. Thus, 176,000 simulation experiments were run in total (**Table 1**). In 175,877 of those simulations, the duration of the infection did not exceed the observation period and were used for further analysis.

#### Analysis

We wanted to determine the probability that a disease outbreak occurs for a certain level of α. The *outbreak probability* was estimated as the proportion of the 2,000 simulation runs, in which the disease spread beyond the starting node. In those cases where the disease spread beyond the starting node, the *outbreak size* was calculated as the total number of infected nodes. In addition, the *outbreak duration* was defined as the number of weeks in which infected nodes occurred. The distribution of the latter two measures was skewed to the right, and thus we give the median and the first and third quartiles (Q1, Q3).

All analyses were conducted using the open-source software *R* version 3.2.1 (41). The package *igraph* (42) was used to generate and analyze the network.

#### RESULTS

#### Descriptors of *G*\*

For this static representation *G*\*, we found an average shortest length of 6.33; the path length between the two most distant nodes (diameter) was 17. The median in-degree, measuring the number of trade partners delivering animals to a certain node was only one, while the median for the number of trade partners a certain holding delivers to (out-degree) was two (**Table 2**). The values for the median ingoing and outgoing closeness centrality were rather similar (**Table 2**), indicating that the number of steps required to reach a certain node equals the number of steps required to reach

any other node from a certain node. The number of shortest paths going through a certain node (betweenness centrality, **Table 2**) showed a high variation, ranging from 0 to 15,166,160.

#### Outbreak Probability

We observed that the outbreak probability is finite, independent of the particular values of transmission probability *p* and network activity α (**Figure 3**). Even for the smallest considered activity values (α = 0.01), the outbreak probabilities was in the region of 5% for all considered transmission probabilities.

TABLE 2 | Minimum, 25% quartile, median, 75% quartile and maximum of the calculated centrality parameters for *G*\*, the static representation of the observed network.

We now focus on the outbreak probability for a transmission probability of *p* = 1, i.e., the worst-case scenario, in which transports of any size spread the infection. In this scenario, a monotonous increase of the outbreak probability with increasing activity was observed. The outbreak probability saturated for larger values of α. More precisely, the outbreak probability was greater than 99% for all α > 0.80. For small and intermediate values of α, it can be observed that even relatively small changes in α had a strong effect on the outbreak probability. Our observed network (α = 0.11) lies in this region. Consequently, small changes in the real system would result in large changes in the outbreak probability.

We now focus on transmission probabilities of *p* < 1. For all considered *p* < 1, a qualitatively similar behavior could be observed. Contrary to the worst-case scenario (*p* = 1), the outbreak probabilities for *p* < 1 did not increase monotonously, but rather showed a maximum. The location of these maxima was shifted to the right for increasing values of *p*. It should be noted that the location of these maxima was relatively close to the activity of the observed network \*.

#### Final Outbreak Size

We now consider the cases where the infection spread beyond the starting node and the corresponding outbreak sizes for different values for α and *p* (**Figure 4**). For the worst-case scenario *p* = 1, the outbreak size increased monotonously with increasing α. The

orange line represents α for the observed pig trade network \*.

possibility that all nodes in the network became infected was only found at this scenario (*p* = 1), but only for very high network activities. For the observed network (α= 0.11), ~15% of the nodes would become infected in the worst-case scenario.

For smaller transmission probabilities (*p* < 1), we observed that outbreak sizes are significantly smaller than in the worst-case scenario. In contrast to the worst-case scenario, the outbreak sizes showed a maximum at approximately α = 0.3.

The authors would like to stress the fact that the outbreak size showed a critical threshold regarding the network activity. This means that there was a *critical activity* αcrit, such that finite outbreaks occurred only if α>αcrit. To estimate αcrit, we calculated the central point between the last value of α below and the first value above the threshold. For transmission probability *p* = 1, we found αcrit = 0.1, and for transmission probabilities of 0.75, 0.5, 0.25, we found αcrit = 0.15. Interestingly, the activity of the observed network was close to the critical region. For *p* = 1, the activity of the observed network \* was only slightly above the critical threshold, whereas for transmission probabilities *p* < 1, the observed network was subcritical. As it is typical for such critical regimes, small changes in the activity result in large changes in the outbreak size (**Figure 4**).

#### Outbreak Duration

Although the shapes of the outbreak durations were similar for different transmission probabilities, we found that the outbreak duration increased with higher transmission probabilities (**Figure 5**). However, for all transmission probabilities, a maximum in the outbreak probability at approximately α = 0.2 could be found, with the exception of *p*= 0.25, where the maximum was at approximately α= 0.3. The reason for these maxima is the existence of two dueling effects. (i) For small α, the outbreak duration correlates with the outbreak size. Outbreaks were typically small here, and increasing α increased the possible number of paths to other nodes. Topological and temporal shortcuts played a minor role here. (ii) For large values of α, the network was likely to form a number of shortcuts, accelerating the spread of a disease.

#### DISCUSSION

In this study, we investigated how the spreading of hypothetical infectious diseases through a trade network is influenced by the networks activity level. For the observed German pig trade network α = 0.11, thus each existing trade connection between two farms was on average active at about 10% of the time during the observation period (using weekly time steps). At this observed low network activity, the chances for a disease to spread beyond the starting node were relatively low, especially for low transmission probabilities (e.g., 10% at *p* = 0.25). Even in the case that an infection spread, the total number of infected farms was for all but the worst-case scenario only about 0.2% of the nodes within the network. Previously, the size of the largest connected component

has often been used as an estimate for the potential final size of an epidemic spreading through a network (43). However, even at the applied worst-case scenario (*p* = 1), the size of the epidemic in our simulation was only a fraction (around 16%) of the total number of nodes for the observed trade data. Using the LSCC would therefore have considerably overestimated the final size of the epidemic. Thus, our results indicate that, at the observed level of the network activity, the threat of large epidemics spreading through the German pig trade network is relatively low, especially for diseases with low transmission rates. However, as we focused in our study on the spread of diseases through the trade network, the actual number of infected farms could be higher due to additional spreading *via* other infection routes (2).

For our analysis, we limited the trade network to the LSCC in order to avoid bias in the results caused by differences in the reachability of the nodes. The observed network activity of the untailored network (α = 0.105) was very similar to the activity in the LSCC (α = 0.112), indicating that changes from the observed network activity would have the same effect on the outbreak probability, -size, and -duration. However, due to the differences in the reachability of the nodes, a much higher variation in the results is to be expected (20).

Interestingly, the German pig trade network seems to be at a rather unstable state, as even small changes in the networks activity level would have a large impact on the spreading of diseases. The main factor that could change the network activity of the German pig trade network is likely to be the farm size. In the last years, the pig production in Germany and other EU countries increased, resulting in larger farm sizes and increased number of traded pigs (1, 44). This would also result in increasing animal transports, which could by archived either by increasing the animals per transport (i.e., edge weights) or by a higher frequency of transports (i.e., higher activity level), whereby the latter would likely have a higher impact on disease dynamics. If an increase in the network activity is to be expected in the long term, the probability for an outbreak and outbreak size are likely to increase, as shown in this study. Considering all three investigated measurements (outbreak probability, final outbreak size, and duration of epidemic), it becomes apparent that an increase in the network activity should be avoided. Further, in order to confine disease spreads, a decrease in the activity of the German pig trade network would be conducive, even if this reduction would only be minor. In our model, a decrease of the activity is realized by random deletion of edges. We assume that a targeted deletion of edges might even have a larger effect (45). From a practical point of view, a reduction in the network activity would mean that animal transports from one farm to another would have to be concentrated to fewer occasions. This also implies that a matching pig production schedule would be necessary, favoring "all-in-all-out" production systems.

The final outbreak size for different network activities shows, as depicted in **Figure 4**, strong similarities to the threshold behavior known from epidemic SIR-type models (46). This epidemic threshold describes a condition above which an epidemic becomes global, while below this threshold only a limited number of nodes become infected (46, 47). To estimate the epidemic threshold in a given network is thus important as it allows predicting the possibility that an infection spreads on a large scale. Hence, it is essential for the planning control and intervention strategies. Different methods exist to identify the epidemic thresholds, with the performance of those methods depending on the topology of the network (48, 49). The results of our study show not only the existence of a threshold but also that its position varies with the transmission probability. Read and colleagues (22) demonstrated for a small-scale human contact network that the encounter rate had a strong effect on the outbreak size at high transmission rates but could find no significant effect at low transmission rates. This concurs with our results, where the effects of the network activity on the outbreak size were most produced at high transmission rates. Again, it seems that the actual activity of the investigated system is close to this threshold value, as even a small increase in the activity level has a large impact on the outbreak size of an epidemic.

The outbreak probability peaked in a region below this threshold for a global epidemic. As the total number of transported animals was kept constant for all network activities, the batch sizes per transport increased, while the frequency of transports decreased at low network activities. Thus, as the edge infection probability *pe* depends on the batch size, the chance to transmit a disease beyond the starting node is rather high at low network activity levels, given the case that a transport occurs. As the number of transports is low at low levels of α, the epidemics are restricted to only a few livestock holdings. On the other hand, a decrease in the observed outbreak probabilities for large values of α can be observed, which can be explained by the fact that the batch sizes are small in this regime. Thus, disease spread was dominated by strong fluctuations in edge infection probability *pe*. These described effects only apply to transmission probabilities <1, as in the case of *p* = 1, the spreading of diseases is independent of the batch size.

As the necessary information was not available to us, we had to made several simplifications for our analysis. Especially, the farm type has already been shown to be an important factor in the spreading of the disease in animal trade networks (15). The farm type defines how long animals remain at a certain node. It is likely, that if the network activity would change due to an overall increase or decrease in the German pig production, the change in the activity of the individual trade connections would be irregular and vary according to the type of the source and the receiving node. This would be an important point to consider in further studies, as heterogeneous waiting times have been shown to influence the spread of diseases in networks (50, 51). For our simulation, we neglected within-herd transmission dynamics as well. Within-herd transmission depends not only on the specifics of a disease but is also influenced by several external factors that were not available to us (e.g., farm size or biosecurity measures on the farm level). The numbers of infected animals within a farm vary over time (52), and it is unlikely that all animals are simultaneously infected over a certain time period, as assumed in our simulation. Consequently, the presented results could overestimate the probability of a disease outbreak and the size of the epidemic. For our model, we assumed that the epidemiological status of the farms does not alter the trade contact structure. This applies to rather harmless diseases, like porcine reproductive and respiratory syndrome (PRRS), porcine circovirus type 2 (PCV2), or MRSA. However, depending on the severity of a disease, trade connections with an infected farm could cease. The withdrawal of trade connections would not be instantaneous but depend on various factors like incubation period or the occurring of clinical symptoms, resulting in a high variation between the time of infection of a farm and the potential termination of trade connections. Thus, the more likely a disease results in trade restrictions and the faster those restrictions are applied, the more our model is prone to overestimate the size of an epidemic. In case of an outbreak of a severe disease, trade connections could change due to the targeted implementation of trade restrictions by veterinary authorities. However, the extent of trade restrictions often differs between countries. For instance, during the bluetongue virus outbreak in Europe starting in 2006, trade restrictions in France were directed to specific areas (53), while in Germany, as well as in Austria and Swiss, the whole country was declared a single restriction zone at an early stage of the epidemic (54–56). Thus, if the whole country is declared a single restriction zone, the within-country trade network would likely show only marginal changes. The effect of lowering the contact rate on outbreak probability, -size, and -duration is shown in this analysis, but the implementation of trade restrictions directed to specific areas could lead to different dynamics.

In our study, we presented the *network activity* as a new indicator value for networks. With this parameter, it is possible to investigate how changes in the mean frequency in the activation of existing trade connections can affect the spread of diseases. By setting the total trade volume constant, as we did in this study, it was possible to differentiate between effects of trade frequency and trade volume. There are two specific characteristics of α: first, it is designed to be a characteristic of a temporal network. It has been shown that several network parameters drawn from a static network correlate with standard epidemiological parameters. Especially in networks with a right-skewed degree distribution, as we found for the pig trade network, nodes with a high degree can play an important role in the spreading of diseases (8, 14, 16). However, the frequency of trade links cannot be represented by a static network; static networks generated from different levels of α would be identical and thus network measurements (like centrality measurements) would be identical as well. As static networks do not take the temporal causality of the paths into account, results drawn from such static representations can be problematic. For example, it has been shown that compared to a temporal network, its static representation overestimates the size of a disease outbreak (20). Thus, in the last years, measurements for temporal networks have been developed (20, 57), and their relation to disease spread, however, remains to be investigated. In comparison with most of those measurements, the calculation of the *network activity* is simple, as it is obtained from the total number of edges in the static and the temporal network. Second, the *network activity* is a measurement for the state of the whole network and not for single nodes. It can be used as a measurement of how well a temporal network is described by its static representation. An α = 1 would be equal to network, where each existing trade link is active at all time steps, thus the static representation would be true at any time. The more closely the network activity is to one, the more accurate is its static representation. Still, for now, we would like to suggest carefulness in applying the results to other networks. While the general pattern is likely to stay the same, the exact location of the maxima/threshold of the investigated parameters could vary. Further, when comparing the *network activity* of different networks, care must be taken to use the same time period and time steps, as α changes with those two values.

In this study, we could demonstrate that the network activity α is an important factor in evaluating the effects of a disease spread in the German pig trade network. We would like to propose applying this indicator number to other networks used to demonstrate the spread of disease or other malicious agents as well, as the networks' activity is likely to have a strong impact on the spreading.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

This work was designed by TS, KL, and HL. BP, HL, and KL processed the raw data, and KL performed the data analysis. The results of the analysis were interpreted by KL, HL, and TS. KL, BP, and HL drafted and wrote the manuscript. All authors revised the manuscript and approved to the final version.

#### ACKNOWLEDGMENTS

We thank A. Fröhlich and C. Firth for their comments, which substantially helped to improve the manuscript.

#### FUNDING

This work was supported by Federal Ministry of Education and Research (BMBF, Bonn, Germany) research grant 13N11208 (KL, HL) and by the project VET-Austria (BP), a cooperation between the Austrian Federal Ministry of Health, the Austrian Agency for Health and Food Safety, and the University of Veterinary Medicine Vienna.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2016 Lebl, Lentz, Pinior and Selhorst. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Prediction of Pig Trade Movements in Different European Production Systems Using Exponential Random Graph Models

*Anne Relun1,2 , Vladimir Grosbois2 , Tsviatko Alexandrov3 , Jose M. Sánchez-Vizcaíno4 , Agnes Waret-Szkuta5 , Sophie Molia2 , Eric Marcel Charles Etter2 and Beatriz Martínez-López1 \**

*1Center for Animal Disease Modeling and Surveillance (CADMS), VM: Medicine and Epidemiology, University of California Davis, Davis, CA, USA, 2CIRAD, UPR AGIRs, Montpellier, France, 3Animal Health and Welfare Directorate, Bulgarian Food Safety Agency, Sofia, Bulgaria, 4Animal Health Center (VISAVET), Animal Health Department, Veterinary School, Complutense University of Madrid, Madrid, Spain, 5 INRA, INP, ENVT, UMR 1225, IHAP, Université de Toulouse, Toulouse, France*

#### *Edited by:*

*Tariq Halasa, Technical University of Denmark, Denmark*

#### *Reviewed by:*

*Anette Boklund, Technical University of Denmark, Denmark Gaëlle Nicolas, Université libre de Bruxelles, Belgium*

#### *\*Correspondence:*

*Beatriz Martínez-López<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Veterinary Epidemiology and Economics, a section of the journal Frontiers in Veterinary Science*

*Received: 06 October 2016 Accepted: 15 February 2017 Published: 03 March 2017*

#### *Citation:*

*Relun A, Grosbois V, Alexandrov T, Sánchez-Vizcaíno JM, Waret-Szkuta A, Molia S, Etter EM and Martínez-López B (2017) Prediction of Pig Trade Movements in Different European Production Systems Using Exponential Random Graph Models. Front. Vet. Sci. 4:27. doi: 10.3389/fvets.2017.00027*

In most European countries, data regarding movements of live animals are routinely collected and can greatly aid predictive epidemic modeling. However, the use of complete movements' dataset to conduct policy-relevant predictions has been so far limited by the massive amount of data that have to be processed (e.g., in intensive commercial systems) or the restricted availability of timely and updated records on animal movements (e.g., in areas where small-scale or extensive production is predominant). The aim of this study was to use exponential random graph models (ERGMs) to reproduce, understand, and predict pig trade networks in different European production systems. Three trade networks were built by aggregating movements of pig batches among premises (farms and trade operators) over 2011 in Bulgaria, Extremadura (Spain), and Côtes-d'Armor (France), where small-scale, extensive, and intensive pig production are predominant, respectively. Three ERGMs were fitted to each network with various demographic and geographic attributes of the nodes as well as six internal network configurations. Several statistical and graphical diagnostic methods were applied to assess the goodness of fit of the models. For all systems, both exogenous (attribute-based) and endogenous (network-based) processes appeared to govern the structure of pig trade network, and neither alone were capable of capturing all aspects of the network structure. Geographic mixing patterns strongly structured pig trade organization in the small-scale production system, whereas belonging to the same company or keeping pigs in the same housing system appeared to be key drivers of pig trade, in intensive and extensive production systems, respectively. Heterogeneous mixing between types of production also explained a part of network structure, whichever production system considered. Limited information is thus needed to capture most of the global structure of pig trade networks. Such findings will be useful to simplify trade networks analysis and better inform European policy makers on risk-based and more cost-effective prevention and control against swine diseases such as African swine fever, classical swine fever, or porcine reproductive and respiratory syndrome.

Keywords: ERGM, network modeling, livestock contact networks, risk-based surveillance, infectious diseases

#### INTRODUCTION

Movements of animals play a key role in the spread of several major infectious diseases, like foot-and-mouth disease, classical swine fever, or African swine fever (1–3). Therefore, detailed data on livestock movements may help to better simulate transmission dynamics and identify areas, periods, and farms that are more likely to spread the diseases and could be targeted to improve surveillance and control strategies (4, 5). However, one of the challenges of using livestock movement data to support decision-making in preventive veterinary medicine is the limited availability of timely and updated records on animal movements and, if available, the massive amount of data that have to be processed. This is particularly challenging when considering diverse and, sometimes, epidemiologically complex, production systems, such as backyard or extensive systems, where the information may not be frequently collected and accessible. Models of livestock movement networks based on holding characteristics and past-temporal observed networks could be useful to simplify real-world networks and to predict disease spread even in backyard or extensive environments.

Pig trade movements can be represented as a network, consisting of a set of nodes (here the pig premises) connected by links (also called edges) representing movements of pigs between them. These networks are not strictly identical from 1 year to the following, but their structural properties, which impact disease dynamics, are likely to be stable over time (6, 7).  For example, some premises may be more likely to trade with each other due to geographical proximity or because they belong to the same pig company [selective mixing or homophily, see Morris (8)]. Some particular types of premises may also be more likely to trade with a high number of premises (attributes that influence degree). Finally, if a trading partner B of premises A trades with a third premises C, this might encourage A to trade with C (structural balance effect).

The first statistical models developed to evaluate which processes lead to observed network structures were quite simple. They only addressed relational reciprocation [i.e., mutuality; see Holland and Leinhardt (9)] or assortative mixing (8). The recent developments of exponential random graph models (ERGMs), also known as *p\** models (10), offer possibilities to better capture the complexity of real-life networks (11). This family of models assumes that the observed network is only one realization among many potential networks with similar characteristics and that the probability that a link exists is a logit-linear function of predictors that reflect node characteristics, link characteristics, and network structural properties (10, 12, 13). Although they were developed to handle the inherent non-independence of network data, the results of ERGMs are interpreted in similar ways to logistic regression, making this a very useful method for examining contact networks in the context of epidemiology.

The aims of this paper were to use ERGMs to (1) develop models that reproduce observed pig trade networks; (2) understand the mechanisms that underlie the organization of pig trade networks; and (3) predict trade networks structures in three different European pig production systems (i.e., industrial, extensive, and backyard). Results of this study are intended to inform the design of prevention and control programs for swine diseases such as African swine fever, classical swine fever, or porcine reproductive and respiratory syndrome under diverse epidemiological scenarios and pig productions systems in Europe.

#### MATERIALS AND METHODS

#### Data Collection and Network Construction

Three areas were selected to represent different European pig production systems: Bulgaria, where most premises raise pigs for own consumption; the autonomous community of "Extremadura," which is the cradle of extensive Iberian pig production in Spain; and the department of "Côtes-d'Armor," which is the French department with the highest concentration of industrial pig premises.

Data on pig movements and premises characteristics were obtained from national databases, through Bulgarian Food Safety Agency in Bulgaria, the professional database of swine (La Base de Données Professionnelle Porcine—BDPORC) in France, and the Ministry of Agriculture, Food and Environment (MAGRAMA) in Spain. The year 2011, which was common in all databases, was retained for analysis. The premises characteristics available were the classification or type of farm (described in the next two sentences), the size of premises (i.e., number of sows, weaners, and finishers on farm), the type of housing system (i.e., indoor or outdoor), the geographical coordinates, and the pig company number (only for France). In Bulgaria, pig farms were classified as small producers (<10 pigs kept for own consumption), type B (medium-size: 10–500 pigs; with low biosecurity level: access to other pigs or feral pigs, use of swill feeding, no fences around the holdings, and/or no disinfection at the entrance and exit of buildings), type A (medium-size, high biosecurity level), or industrial farms (large size: >500 pigs; high biosecurity level) (14, 15). Traditionally and outdoor-raised East Balkan pig herds are also found in the South East of Bulgaria. For Spain and France, pig farms were classified as multipliers (premises that produce breeding stocks and semen), farrowing farms, farrow-to-finish farms, finishing farms, or small producers. Small producers for Spain were defined as those that produce pigs for own consumption, whereas for France were those with ≤4 pigs. Traders, collection centers, markets, fairs, and stopping points (i.e., or staging point: locations used to feed, water, rest, accommodate, care for, and dispatch animals in transit before arriving to their final destinations) were considered as trade operators. Because of the dead-end characteristics of slaughterhouses, these premises were excluded from analysis.

For each area, yearly networks (i.e., using year 2011) were built, the nodes being all pig premises of the study areas, even those that were not trading pigs during the study period. Movement data were aggregated over the study period, and a direct link was drawn whenever a shipment of pigs occurred between the corresponding premises. Movement imported from or exported to outside areas was excluded from the analyses.

#### The ERGMs

Exponential random graph models specify the probability of any random network **Y** given a set of *n* nodes and their attributes as in Eq. 1.

$$P\_{\theta} \left( \mathbf{Y} = \boldsymbol{\mathcal{Y}} \mid n \text{nodes} \right) = \left( \frac{1}{c} \right) \exp \left( \sum\_{k=1}^{K} \theta\_{k} z\_{k} (\boldsymbol{\mathcal{Y}}) \right) \tag{1}$$

The *zk*(*y*) terms represent model covariates, any set of *K* network statistics calculated on the *y* observed network and hypothesized to affect the probability of this network forming. The model covariates can include network parameters that account for the frequency of occurrence of certain network configurations (e.g., two-path, triangles), as well as node or edgewise covariates like the pig company to which a premise belongs or the distance between two premises, respectively. The θ coefficients estimate the strength of the effect of each covariate. The denominator *c* represents the normalizing constant, which correspond to the sum of exp ( θ ) *k k k z y* =1 K ( ) ∑ over all possible networks with *n* nodes.

Because ERGMs' calculation time dramatically increased with the increase of network size, it was decided to exclude isolates, i.e., pig premises that did not trade with other premises, from the small-scale productions system (Bulgaria, initially 28,729 premises, of which 95.3% were isolated premises).

#### Model Specification

First, an exploration of network data was undertaken, with the computation of several local topological measures (number of isolates, triangles, degree distribution, etc.) and of mixing matrices for premises' attributes (16). Specifically, we computed the number of nodes, the network density, the percentage of isolates, the clustering coefficient, and the mean and range of in-degree and out-degree centrality measures [e.g., Ref. (5)]. Network graphs were plotted, with the nodes colored according to nodes' attributes, to better visualize the selective mixing patterns.

Based on this exploration, several network statistics were chosen to represent hypothetic rules for trade movements (**Table 1**). *L*(*y*) captures the density of the observed network *y*. *Mi,v,a*(*y*), *Mo,v,a*(*y*), *Ha,v*(*y*), *Ua,v*(*y*)*, Sa,v*(*y*), and *E*(*y*) are attributespecific terms that capture the way in which premise attributes structure trading patterns, where *a* represents the attribute (e.g., housing system) and *v* the level (e.g., indoor, outdoor). The main effects, *Mi,v,a*(*y*) and *Mo,v,a*(*y*)*,* allow variation in the propensity of a premise to form incoming and outgoing edges according to the level of an attribute characterizing this premise. *Ha,v*(*y*) models a tendency of edges to occur between premises belonging to the same attribute level that varies among attribute levels (hereafter referred to as *differential homophily*), while *Ua,v*(*y*) models a uniform tendency of edges to occur between premises belonging to the same attribute level (hereafter referred to as *uniform homophily*). *Sa,v*(*y*) accounts for variation in the occurrence of edges according to the levels of an attribute characterizing each of two premises (hereafter referred to as *selective mixing*). *E*(*y*) captures variation in the propensity of premises to form links according to the Euclidean distance in km to other premises.

Table 1 | Network statistics used to fit the exponential random graph models of pig trade networks.

*a Some statistics use attribute-specific terms where a and v represent the attribute and level, respectively. The observed network is represented by y and the scale parameter by* α*.*

*A*(*y*) and *I*(*y*) model the tendency of premises to form unidirectional links or no links, respectively. The terms *gwdsp*(*y,* α), *gwesp*(*y,* α), *gwid*(*y,* α), and *gwod*(*y,* α) are related to local structures and represent the parametric forms of the alternating twopaths, clustering (alternating *k*-triangles) and in- and out-degree distributions, respectively. A fixed value of 0.5 was adopted for the scale parameter α in these terms (11).

The Markov Chain Monte Carlo (MCMC) algorithm was used to estimate the maximum likelihood for the θ coefficients included in models (12). The MCMC chain is intended to step around the sample space of possible networks, selecting a network at regular intervals to evaluate the statistics in the model. For each MCMC step, *n* (*n* = 1 in the simple case) toggles are proposed to change the dyad(s) to the opposite value. A chain burn-in of 105 toggles, an MCMC sample size of 104 , and an interval between successive samples of 103 was fixed for these models.

#### Model Selection and Goodness of Fit

For each area, four models were built: (1) a simple Bernoulli model that only includes the number of edges; (2) a model with edges and statistics based on nodal attributes (hereafter called "edge + attribute" model); (3) a model with edges and structurerelated statistics ("edge + network statistics" model); and (4) a model with edges, nodal attributes, and structure-related statistics ("edge + attributes + network statistics" model).

For the "edge + attribute" and "edge + network statistics" models, univariable analyses were performed first. The terms (i.e., attributes and network statistics) were then added one by one, until the best model fit was obtained. The fourth model was based on the final "edge + attribute" model, and network statistics terms were added one by one manually until the best model fit was obtained.

Three approaches were used to examine goodness of fit of the models: (1) check for model convergence and degeneracy; (2) comparison of Akaike information criteria; and (3) comparison of goodness of fit plotting for higher order statistics (11). For this purpose, four sets of statistics were used: the in- and out-degree distributions, the geodesic distance distribution, and the edgewise shared partner distribution, which reflects the clustering of the network (17). These statistics were chosen because of their impact on disease spread dynamics (18). Finally, plots of simulated networks were visually compared to the plot of the observed networks.

All analyses were conducted in R (19) using the "statnet" suite of packages (20, 21).

#### RESULTS

A total of 7,811 out of the 45,224 premises keeping pigs (i.e., 17.3%) were actively moving pigs during 2011. Description of the pig industry demographics (i.e., number of premises for each type of farm), pig trade (to/from different types of farm), and topological characteristics in Côtes-d'Armor (France), Bulgaria, and Extremadura (Spain) in 2011 are presented in **Tables 2**–**4**.

The inclusion of both nodal attributes and network configurations statistics provided the best fit to the data (**Tables 5**–**7**; **Figures 1** and **2**). Selective mixing between premises according to their type of production appeared to be an important mechanism of pig organization, whichever system considered (**Tables 5**–**7**). In addition to this mechanism, the other mechanisms related to premises characteristics that impacted the most on trade organization were belonging to the same pig company, the tendency of outdoors premise to trade with outdoor premises, and the geographical location of pig premises, in the intensive, extensive, and small-scale production systems, respectively (**Tables 5**–**7**). Network statistics on dyadwise and edgewise shared partner distributions, as well as on in- or out-degree distributions were

*a For Bulgaria only: industrial farm (large size:* >*500 pigs, high biosecurity level farm); type A farm (medium-size: 10–500 pigs, high biosecurity level); type B farm (medium-size: 10–500 pigs, low biosecurity level); east Balkan pigs (traditional outdoor pig herds). Small producer for Bulgaria:* <*10 pigs kept for own consumption. bSmall producers were defined for Spain as those farms (any size) that produce pigs for own consumption; for France were those with* ≤*4 pigs and for Bulgaria were those with* <*10 pigs kept for own consumption. NA, not applicable/not available.*

*IQR, interquartile range.*

*a Premises that sent or received pigs in 2011.*

#### Table 4 | Topological statistics of the pig trade networks in 2011.

*kin, in-degree; kout, out-degree.*

*a Isolates were excluded for Bulgaria: initially there were 28,729 premises, of which 95.3% were isolated premises.*

Table 5 | Parameter coefficients and fit for the four exponential random graph models (ERGMs) of pig trade in a small-scale production system (Bulgaria).

*a \*\*\** <*0.001; \*\** <*0.01; \** <*0.05.*

*bE, east; NW, north-west; S, south; SW, south-west.*

*c SP, small producer; IN, industrial; TA, type A; TB, type B.*

needed to fit the models (**Tables 5**–**7**). These statistics better reflected the clustering of the observed networks and allowed to reproduce well the observed global network properties (**Figures 1** and **2**). This can be clearly observed in the goodness of fit diagnostics plot (**Figure 2**), where the value of all the observed network statistics (solid black line) is only well captured by the distribution of values of the simulated networks (underlying boxplots) generated with the final ERGM model (i.e., model with edges + attributes + network statistics).

#### DISCUSSION

Exponential random graph models were used to represent, understand, and predict pig trade networks structures from different European production systems, with predominantly small-scale, extensive, or intensive pig producers. Such information improved our understanding of the processes that govern the organization of pig trade and can further be used to better inform European policy makers on prevention and control

Table 6 | Parameter coefficients and fit for the four exponential random graph models (ERGMs) of pig trade in an extensive production system (Spain—Extremadura).

*a \*\*\** <*0.001; \*\** <*0.01; \** <*0.05.*

*bIn, indoor; Out, outdoor.*

*c U, multipliers; FA, farrow farms; FF, farrow-to-finish farms; FI, finishers; SP, small producers.*

measures against swine diseases such as African swine fever, classical swine fever, or porcine reproductive and respiratory syndrome.

Rapidity of targeted action during the initial phase of an outbreak is fundamental to effectively curtail the transmission and minimize the disease burden. At this time, movements of animals have not been banned, and it is thus relevant to use "peace-time" movement networks to compare different control strategies. Until recently, variability in contact patterns was mostly approached in epidemic models by combining probabilities of contact between premises according to their type of production and to the distance between premises (22, 23). These efforts may fail to capture the structural properties of livestock trade networks that will impact diseases dynamic as well as their spatial spread (24–26). The use of ERGMs to model pig trade networks allows capturing both network topology and complex behaviors that depend on various premises characteristics. They may thus help to generate more realistic networks that may be used to study diseases spread, identify premises that could be targeted for risk-based surveillance, early detection, and rapid control of diseases, and compare different control strategies (27–29). Indeed, by simulating diseases spread on several simulated networks, we could identify some farms that are frequently and early infected and thus that should be targeted to provide timely and accurate indications of epidemic activity (30, 31). These networks could also be simulated to

Table 7 | Parameter coefficients and fit for the four exponential random graph models (ERGMs) of pig trade in an intensive production system (France—Côtes-d'Armor).

*a \*\*\** <*0.001; \*\** <*0.01; \** <*0.05.*

*bFor readability, not all values for selective for the pig companies are shown; NC, no company.*

*c MU, multipliers; FA, farrow farms; FF, farrow-to-finish farms; FI, finishers; SP, small producers; TR, trade operators.*

assess the effectiveness of compartmentalization or zoning, strategies that might be efficient to prevent disease's spread without disrupting pig trade (32, 33).

Exponential random graph models developed in this study also improved our understanding of the drivers of pig trade in different production systems. Geographic mixing patterns

strongly structured pig trade organization in the small-scale production system, whereas belonging to the same company, or keeping pigs in the same housing system appeared to be key drivers of pig trade, in intensive and extensive production systems, respectively. As expected, the specialization and organization of pig production also explained a part of trading

behaviors, illustrated by the heterogeneous mixing between types of production. This mechanism was however less important than the geographical location of premises and as important as belonging to the same pig company in the small-scale and intensive production systems, respectively. Geographical proximity did not appear to play a role in the intensive system, whereas it was significant in the small-scale producer system. Unfortunately, model degeneracy occurred when trying to include the distance between premises as a covariate in the extensive production system (Spain—Extremadura), preventing a conclusion to be drawn on the impact of this covariate in this production system.

Finally, this study revealed that the inclusion of nodal attributes was necessary to represent the mixing patterns, but it was not sufficient to reproduce the great clustering observed in the pig trade networks, which could only be represented when adding additional statistics on local network configurations. These statistics reveal some features of these networks, such as the propensity of trade to have a short path length (negative coefficient for the degree distribution terms). Some social behavior, economic factors, or unobserved covariates, which may differ between countries, may also have driven the choices of farmers for trading partners (e.g., pig prices, road distribution, traditions, or cultural practices, etc.). This may explain the increased clustering as represented by the positive coefficient for the *gwesp* term.

Until recently, problems of degeneracy and computational intractability for large network sizes limited the use of ERGMs in epidemiological modeling (11, 34). Indeed, ERGMs have been mainly used on small networks to understand the factors driving human social behaviors (35, 36) and have sometimes been applied on disease transmission modeling (37). Ortiz-Pelaez et al. (38) were the first to introduce ERGMs in preventive veterinary medicine, using this method to understand the factors driving livestock trade in a small network of villages in Ethiopia. In the present study, the use of new parameters that limit degeneracy problems (12) allowed us to obtain statistical models with a good fit to the large-size observed networks. Isolates always depends on the spatial and temporal "frontiers" that are decided and on the exchange with farms outside these "frontiers." For Extremadura, most of movements were inside the region (i.e., from the 9,544 isolates, 174 of them sent pigs to farms outside Extremadura and 73 received pigs from farms outside Extremadura); therefore, most of the isolates (97.4%) can be considered as true isolates during the study period (i.e., not movements within the region and not trading with other regions). For Côtes-d'Armor, there was a lot of exchange with other departments, and only 47% of the 489 isolates can be considered as true isolates (i.e., 57 sent pigs to farms outside Côtes-d'Armor and 202 received pigs from farms outside the region). Therefore, the scale to simulate networks for disease spread models should consider areas where almost all movements are inside these areas. For Bulgaria, the entire country was evaluated, and all were considered to be true isolates; however, it was not feasible to include these isolates in the ERGM due to memory limits, and therefore the model produced here might not fully represent the true pig movement network in Bulgaria. Further studies should be conducted to validate this network once computational difficulties to fit ERGM's to large networks are solved.

Since implementation of Regulation (EC) no 1760/2000 of the European parliament, recording of livestock movements between premises is mandatory, making data on pig trade movements available, at least in the main producing countries in the EU. However, there are no standards on the definition of different types of premises (e.g., backyard) or on other premises attributes. The scales of the networks considered in this study were also different, being the national level for Bulgaria and the regional level for France and Spain, to better study very specific production systems. Therefore, though the analyses and results in the different settings intended to illustrate the applicability and usefulness of the approach in the predominant swine production systems in the EU, mechanisms and rules that govern trade organization in the different study populations are not fully comparable.

Several studies showed that, in addition to the topology of a contact network, heterogeneity in the weight of edges and temporal network dynamics had a strong influence on diseases spreading (39, 40). Tools to model such networks are still under development, and their application is currently limited by the size of the networks modeled (41–43). In the next few years, these methods could be promising tools to improve our representation of real-world livestock trade networks.

#### CONCLUSION

This study is one of the very first to illustrate the usefulness of ERGMs to understand and simulate livestock trade networks under different European production systems, specifically small-scale, extensive, and intensive swine production systems. Depending on the production system, some premises characteristics, such as their geographical location, type of production, belonging to a pig company or housing system, were key drivers of pig trade, but adding statistics on local network configurations was necessary to accurately capture the great clustering observed in all pig trade networks. These models offer a framework to simulate realistic pig trade networks that may be included in epidemic models to compare different control strategies against major swine diseases such as African swine fever, classical swine fever, or porcine reproductive and respiratory syndrome.

#### AUTHOR CONTRIBUTIONS

AR, BM-L, and VG designed the study and developed the R codes. AR and BM-L gathered, cleaned, and verified the data. TA, AW-S, SM, EE, and JS-V contributed to the interpretation and critical discussion of the nature, characteristics, and structure of the data for the different study regions. AR carried out the analyses and wrote the draft of the manuscript. All authors participated in the interpretation and discussion of the results, read, edit, and approved the final manuscript.

#### ACKNOWLEDGMENTS

The research leading to these results received funding from the European Union, Seventh Framework Programme (FP7/2007- 2013) under grant agreement no 311931 (ASFORCE) and Boehringer Ingelheim Vetmedica, Inc. (gift support). The authors are very grateful to BFSA, BDPORC, and MAGRAMA for providing the data.

#### REFERENCES

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer AB and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

*Copyright © 2017 Relun, Grosbois, Alexandrov, Sánchez-Vizcaíno, Waret-Szkuta, Molia, Etter and Martínez-López. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Early Decision Indicators for Foot-and-Mouth Disease Outbreaks in Non-Endemic Countries

*Michael G. Garner1 , Iain J. East1 \*, Mark A. Stevenson2 , Robert L. Sanson3 , Thomas G. Rawdon4 , Richard A. Bradhurst5 , Sharon E. Roche1 , Pham Van Ha6 and Tom Kompas5*

*1Animal Health Policy Branch, Department of Agriculture and Water Resources, Canberra, ACT, Australia, 2 Faculty of Veterinary and Agricultural Sciences, University of Melbourne, Parkville, VIC, Australia, 3AsureQuality Limited, Palmerston, New Zealand, 4 Investigation and Diagnostic Centre and Response Directorate, Ministry for Primary Industries, Wellington, New Zealand, 5Centre of Excellence for Biosecurity Risk Analysis, University of Melbourne, Parkville, VIC, Australia, 6Crawford School of Public Policy, Australian National University, Acton, ACT, Australia*

#### *Edited by:*

*Tariq Halasa, Technical University of Denmark, Denmark*

#### *Reviewed by:*

*Kimberly VanderWaal, University of Minnesota, USA Preben William Willeberg, Technical University of Denmark, Denmark*

> *\*Correspondence: Iain J. East*

<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Veterinary Epidemiology and Economics, a section of the journal Frontiers in Veterinary Science*

*Received: 26 August 2016 Accepted: 17 November 2016 Published: 30 November 2016*

#### *Citation:*

*Garner MG, East IJ, Stevenson MA, Sanson RL, Rawdon TG, Bradhurst RA, Roche SE, Van Ha P and Kompas T (2016) Early Decision Indicators for Foot-and-Mouth Disease Outbreaks in Non-Endemic Countries. Front. Vet. Sci. 3:109. doi: 10.3389/fvets.2016.00109*

Disease managers face many challenges when deciding on the most effective control strategy to manage an outbreak of foot-and-mouth disease (FMD). Decisions have to be made under conditions of uncertainty and where the situation is continually evolving. In addition, resources for control are often limited. A modeling study was carried out to identify characteristics measurable during the early phase of a FMD outbreak that might be useful as predictors of the total number of infected places, outbreak duration, and the total area under control (AUC). The study involved two modeling platforms in two countries (Australia and New Zealand) and encompassed a large number of incursion scenarios. Linear regression, classification and regression tree, and boosted regression tree analyses were used to quantify the predictive value of a set of parameters on three outcome variables of interest: the total number of infected places, outbreak duration, and the total AUC. The number of infected premises (IPs), number of pending culls, AUC, estimated dissemination ratio, and cattle density around the index herd at days 7, 14, and 21 following first detection were associated with each of the outcome variables. Regression models for the size of the AUC had the highest predictive value (*R*2 = 0.51–0.9) followed by the number of IPs (*R*2 = 0.3–0.75) and outbreak duration (*R*2 = 0.28–0.57). Predictability improved at later time points in the outbreak. Predictive regression models using various cut-points at day 14 to define small and large outbreaks had positive predictive values of 0.85–0.98 and negative predictive values of 0.52–0.91, with 79–97% of outbreaks correctly classified. On the strict assumption that each of the simulation models used in this study provide a realistic indication of the spread of FMD in animal populations. Our conclusion is that relatively simple metrics available early in a control program can be used to indicate the likely magnitude of an FMD outbreak under Australian and New Zealand conditions.

Keywords: FMD, early decision indicators, vaccination, simulation models, decision-support, regression analysis

#### INTRODUCTION

Disease managers are faced with a number of challenges when deciding on the most effective disease control strategy to implement in an exotic animal disease outbreak. Foot-and-mouth disease (FMD) is particularly challenging given its wide range of host species, potential for rapid spread, and serious socioeconomic consequences. For countries such as Australia and New Zealand, FMD represents the most serious threat to their livestock industries. A recent study estimated the 2013 value of total direct economic loses over 10 years for a large multi-state outbreak of FMD in Australia at USD 47 billion (1). Animal products constitute a significant proportion of New Zealand exports, and the provisional results of recent modeling of the economic impacts of a large FMD outbreak in New Zealand have estimated net 2014 GDP losses over an 8-year period to be between USD 13 and 17 billion (2). Consequently, Australia and New Zealand invest considerable resources in preparedness and planning for emergency animal disease outbreaks, including maintaining vaccine banks for FMD. Despite recent changes to contingency plans to recognize that vaccination could be an important component of an FMD control program, it is unclear how or when vaccination should be used, and if it is used, how vaccinated animals should be managed once an outbreak has been resolved.

Modeling studies carried out in Australia (3–5) and overseas (6–8) have shown that vaccination is effective in reducing the duration and/or size of FMD outbreaks in situations where disease is widespread, where there is a high rate of spread or the resources for stamping out are limited. Reports suggest that early vaccination may have been beneficial in eradicating the disease earlier than was the case with recent FMD outbreaks in Korea (9) and Japan (10). Thus, vaccination is increasingly being recognized as a potential useful tool to assist in containing and eradicating FMD outbreaks in countries where the disease is not endemic. However, while vaccination may contribute to earlier eradication of the disease, it will be associated with additional costs – keeping vaccinated animals in the population will delay the period until FMD-free status is regained under current World Organization for Animal Health guidelines (11) and add additional complexity to post-outbreak surveillance programs. These issues are of particular concern for countries with significant exports of livestock and livestock products because, under current conditions, the use of vaccination and the presence of FMD vaccinated animals in the population could be expected to cause significant market access difficulties.

From a planning and management perspective, it would be useful to have access to decision support tools that take into account the information that would be available to disease managers early in an outbreak to provide an indication of the potential severity of the outbreak that could ensue. This would enable decisions on specific measures like vaccination to be made at a time when they are likely to be most effective.

McLaws and Ribble (12) documented the relationship between the interval (in days) from incursion to detection and epidemic size [expressed as the total number of infected premises (IPs)] for 24 FMD outbreaks in non-endemic countries that occurred between 1992 and 2003. They did not find a direct relationship between time to detection and total number of IPs or total animals culled for disease control, concluding that the movement of animals through markets was the most critical factor contributing to large outbreaks. Sarandopoulos (13) conducted a review of 125 FMD epidemics in non-endemic temperate countries reported to the OIE between January 1, 2005 and December 31, 2013 to identify associations between epidemic size/duration and early outbreak explanatory variables. The explanatory variables assessed in this study included susceptible animal densities, weather conditions at the time of detection, the number of IPs detected in the first 7 days, and the size of the area under control (AUC) at 7 days (based on a convex hull calculation). In total, ten candidate explanatory variables were tested for their association with epidemic size and duration using a zero-inflated negative binomial regression model. Cattle density, pig density, and the number of IPs at day 7 post-detection were all positively associated with epidemic size while increased average temperature in the month of detection was associated with "smaller" outbreaks.

Using data from the outbreak of FMD that occurred in the UK in 2001, first fortnight incidence (FFI), i.e., the cumulative number of new FMD-IPs found in the first 2 weeks of the response, was found to be a useful predictor of the size and duration of outbreaks at the regional and national scale (14, 15). The larger the number of detected herds within the first 2 weeks, the higher the risk of the large outbreak. Halasa et al. (16) extended the approach of Hutber et al. to incorporate the first fortnight spatial spread (FFS) as well as FFI (which they renamed first fortnight outbreaks – FFO, since a true incidence rate is not actually calculated) – in a simple decision tool using simulated FMD outbreaks. In terms of outcome, in addition to the number of IPs and outbreak duration, they also considered the size of the AUC and costs. Halasa and colleagues found good correlations between FFO and FFS and all of the outcome variables, indicating that both FFO and FFS have the potential as predictors of epidemic outcomes. They also found that the type of index herd was a significant predictor of epidemic outcome.

The combined work of Hutber et al. (15), Halasa et al. (16), and Sarandopoulos (13) indicates that information available early in an outbreak can be used to make inferences about the potential severity of an FMD outbreak and could perhaps be incorporated into decision support tools. However, one of the concerns is that FFO and FFS are quite simple parameters that are likely to be sensitive to outbreak management response, in particular the effectiveness of the surveillance/reporting system. For example, while a low FFO may be indicative of a limited spread and small number of infected places, it could also be indicative of the adequacy of resources to undertake surveillance and tracing. In addition, based on the work of McLaws and Ribble (12) and Sarandopoulos (13), other factors such as animal densities at the location of the index premises and involvement of animal markets may also be important.

With this background, this study was undertaken to identify characteristics measurable during the early phase of a FMD outbreak that might be useful as predictors of the severity of an FMD epidemic (expressed as the total number of infected places, outbreak duration, and the total AUC). The study also aimed to assess how robust findings were across different incursion scenarios and between different production and management systems. A key point is that in this study simulation models of FMD were used to generate a series of outbreaks listing incident infected places over time and geographical space. Regression approaches were then used to identify characteristics measurable during the early phase of a simulated outbreak that might be useful as predictors of the total number of infected places, outbreak duration, and the total AUC predicted by each simulation model. The inferences drawn from this study are dependent on the strict assumption that each of the simulation models used in this study provide a realistic indication of the spread of FMD in animal populations.

#### MATERIALS AND METHODS

A modeling study was undertaken to test a range of explanatory variables as predictors of the total number of infected places, outbreak duration, and the total AUC in an FMD outbreak. The study involved two countries (Australia and New Zealand) and two modeling platforms. Linear regression, classification and regression tree (CART), and boosted regression tree (BRT) analyses were used to assess the association between putative explanatory variables and the three outcome variables.

#### Disease Models

The Australian Animal Disease Spread Model [AADIS (17, 18)] is a hybrid model that simulates the spread and control of FMD in livestock populations at a national scale. AADIS uses the herd as the epidemiological unit of interest and models the spread of disease both within and between herds. Spread of disease within a herd is modeled through a deterministic equation-based model, and between-herd spread is modeled with a spatially explicit stochastic agent-based model. There are five discrete spread pathways in the between-herd model: direct animal movements, local spread (infection of farms within close geographical proximity by unspecified means), indirect contact (*via* contaminated equipment, people, or animal products), animal movements *via* saleyards, and windborne spread.

The model incorporates the attributes and spatial locations of individual farms, saleyards, weather stations, local government areas, and various other features of the environment. For FMD control, AADIS is configured to support the range of mitigation strategies described in Australia's contingency plans for FMD (19) with the effectiveness of these measures dependent on available resources (4).

InterSpread Plus [ISP (20)] is a spatial and stochastic simulation model of infectious disease in domestic animal populations. ISP is a state-transition model meaning that the epidemiological units of interest (farm locations) exist in either the susceptible, infected, or not-at-risk state at any given time. Similar to AADIS, ISP uses a series of user-defined parameters to define the spread of an infectious agent from one farm location to another through local spread, windborne spread, and direct and indirect contacts. Updated movement parameters are informed by findings from recent livestock movement studies in New Zealand (21, 22). Control measures, such as depopulation, vaccination, and movement restrictions, in addition to varying disease surveillance intensity can be simulated, with the ability to carry out each of these activities subject to user-defined resource constraints, similar to the AADIS model.

#### Study Design

Epidemics of FMD in Australia and New Zealand were simulated using the AADIS and ISP models, respectively. A total of 10,000 FMD outbreak simulations were carried out using each model. For each simulation, FMD was introduced into a single livestock farm selected at random within assessed high risk areas for FMD. For Australia, the study area was the whole country, with initial seeding of infection confined to south eastern Australia (**Figure 1A**). South eastern Australia comprises the states of Victoria and Tasmania and parts of New South Wales and South Australia. This area contains a mix of farming enterprises. It is the center of Australia's dairy production and is considered a higher risk area for introduction, establishment, and spread of FMD (23).

The study area for New Zealand comprised the whole of mainland New Zealand, incorporating the North and South Islands. Initial seeding of infection was confined to the Auckland megaregion (Auckland and its three neighboring regional council areas, **Figure 1B**) as it is assumed that the most likely introduction scenario for FMD into New Zealand would involve people or contaminated products seeding infection into livestock in this area. The Auckland mega-region has the largest international air and sea ports. Furthermore, yachts visiting the country are more likely to make landfall in the north.

The following assumptions were used for the Australian and New Zealand FMD models. The time from incursion to first detection was probabilistically determined based on farming systems and expected disease reporting rates in the two countries. For Australia, data on the daily probability of detection and the delay from incursion to first detection were sourced from Martin et al. (24). For New Zealand, data on the daily probability of detection were sourced from Murray and Sanson (25). Outbreak control was based on application of animal movement controls, enhanced surveillance, tracing, and stamping out (i.e., destruction, disposal, and decontamination) on detected IPs. These were applied according to each country's FMD response plan (19, 26). Resources for disease control were based on each country's estimates of expected resources (5). Each model run ended when disease was eradicated or after 1 year, whichever occurred first.

#### Explanatory Variables

Three time points (days 7, 14, and 21 after first detection) were selected, and candidate explanatory variables based on data that would be available to disease managers at these time points were collated: (1) outbreak location: farm and animal densities around the site of first detection; (2) the involvement of markets/ saleyards; (3) measures of the geographic distribution of IPs, as measured by the AUC, and the number of discrete disease clusters; (4) measures of temporal spread, as measured by the number of IPs reported, and the number of traced premises identified; (5) the rate of disease spread, as measured using the estimated dissemination ratio (EDR), calculated using the methods

described by Miller (27) and Morris et al. (28); and (6) adequacy of resources available for control. A description of each of the candidate explanatory variables is provided in **Table 1**.

For each simulated outbreak, we defined three outcome variables: the total number of IPs, outbreak duration (defined as the number of days from first detection until the date on which the last IP was culled), and the total AUC (in km2 ). Modeling results for each country were analyzed separately.

#### Statistical Methods

#### Linear Regression

The Stata/IC statistical package (29) was used for all linear regression analyses. Datasets were imported into Stata and the three outcome variables and each of the explanatory variables checked for normality and log transformed, where necessary, to minimize problems due to non-normality and heteroscedasticity of model residuals (30). Scatterplots of each of the log transformed outcome variables and each of the log transformed candidate explanatory variables were made and the association between variable pairs assessed by superimposing a Lowess-smoothed curve on each plot. After log transformation, all relationships were linear or near linear. Subsequent analyses used two modeling techniques: (1) linear regression modeling using robust estimates to account for non-normally distributed dependent variables (31) and (2) negative binomial regression. It was considered appropriate to use linear regression techniques because the methodology is robust to violations of the requirement for normally distributed dependent variables if the number of observations is large (32, 33). The outputs of the two regression models were similar, so only the results of the linear regressions are presented. The linear regression was preferred because a small proportion of values had excessively large residuals in the negative binomial models.

#### TABLE 1 | Explanatory variables tested.

Candidate explanatory variables were initially tested for unconditional associations with each of the three outcome variables. Explanatory variables that were associated with the outcome variables with *P* < 0.20 were selected for inclusion in the initial multiple regression models. The initial multiple regression model was then reduced step-wise by removing the explanatory variable with the highest *P*Wald value. This process was repeated until all remaining explanatory variables had *P*Wald < 0.05. After the most parsimonious model was developed, all excluded explanatory variables were reassessed by adding them individually back into the model. All biologically plausible, first-order interaction terms were tested, one at a time and retained in the model if the *P*Wald < 0.05 (no interaction terms were retained). The extent of confounding was assessed using the variance inflation factor. No significant confounding was observed in the final models presented.

We acknowledge that the outcome variables measured on a given EDI day (i.e., days 7, 14, and 21 post-detection), which were used as explanatory variables in each model were correlated with their corresponding outcome variable. To investigate this issue further, alternative models were developed where the outcome variable was expressed as (for example) the total number of IPs − IPs identified up to day 14. Using this approach, we identified no substantial differences in the final set of explanatory variables included in the model and the direction and magnitude of the adjusted measure of association between each explanatory variable and the outcome were essentially the same. For this reason, and also to allow our findings to be compared with previous studies (16), we elected to use the total number of infected places, total outbreak duration, and the total AUC as outcome variables in each of the models presented.

While the explanatory variables that remained in each of the Australian and New Zealand regression models differed, those with the most explanatory power (that is, those with the highest beta weight values) were present in both of the country models. For parsimony, a simpler regression model was built using only the explanatory variables that were common to both the Australian and New Zealand models with little or no loss of explanatory power (see Results).

For the linear regression models, the *R*2 value is reported as a measure of the goodness of fit of the model. Based on the regression coefficients estimated for the explanatory variables included in each of the three regression models for each country, predictions of the total number of infected places, outbreak duration, and the total AUC were computed. Several cutpoints (e.g., more or less than 20 IPs) were then arbitrarily selected to divide the model iterations into large and small outbreaks. Two by two contingency tables were constructed to compare the regression model estimates with the actual values of the (classified) outcome variables. These data were then used to calculate negative and positive predictive values for the day 14 model estimates using standard techniques (34).

#### Regression Trees

Acknowledging the possibility of non-linear relationships between the explanatory variables and the three outcome variables used in this study, we used CART and BRT analyses as an alternative approach for identifying associations in these data. CART analysis involves recursively partitioning an outcome variable into two parts based on the value of a given predictor variable that best splits the data.  Predictor variables and their split points are chosen to optimize a given goodness-of-fit criterion, such as minimizing the residual sum of squares (for continuous data). CART analysis is mathematically identical to some multivariable regression techniques, but presents the results in a way that is easily understood by nontechnical audiences.

In contrast to CART, a BRT analysis generates a large number of regression trees based on random samples of the data (35). A BRT model returns a list of predictor variables used to create the splits in each of the trees computed using the randomly sampled data. A relative weight is then calculated for each predictor variable by computing the average number of times the variable was chosen for splitting weighted by the squared improvement to the model from each split and scaled to sum to 100. Larger weights indicate a stronger influence between an explanatory variable and the outcome. The BRT analysis requires the analyst to specify the learning rate and tree complexity. Learning rate controls how much each tree contributes to the model as it develops. In general, smaller learning rates result in better predictions than larger learning rates. Tree complexity sets the number of interactions fitted in the model: a tree complexity of two allows for two-way interactions, three allows for three-way interactions, and so on.

Classification and regression tree analyses were carried out for each of the three outcome variables for the Australian and New Zealand data using the rpart package (36) implemented in R version 3.3.1 (37). The BRT analyses were carried out using the dismo package (38) in R.

#### RESULTS

Of the 10,000 outbreaks that were simulated, FMD did not establish (there was no spread from the seed herd) in 3210 simulations in Australia and 1180 simulations in New Zealand. These simulations were excluded from subsequent analyses. Descriptive statistics of the simulated outbreaks and explanatory variables for Australia and New Zealand are shown in **Tables 2** and **3**, respectively. Descriptive statistics of the outcome variables for the AADIS (Australia) and ISP (New Zealand) models (that is, the total number of infected places, outbreak duration and the total AUC) are shown in **Table 4**.

For the New Zealand (ISP) simulations, the area used for seeding FMD outbreaks had a substantially higher density of cattle (median of 152 head/km2 ) than the areas where FMD was seeded for the Australian (AADIS) simulations (median of 28 head/km2 ).

Compared with the FMD outbreaks simulated by AADIS, ISP simulated relatively high numbers of IPs during the early phase of

*AUC, area under control; EDR, estimated dissemination ratio.*

*a Cattle density: number of cattle/km2 .*

*bSheep density: number of sheep/km2 .*

*c Pig density: number of pigs/km2 .*

*dHuman density: number of humans/km2 .*

*AUC, area under control (km2 ); EDR, estimated dissemination ratio.*

*.*

*a Cattle density: number of cattle/km2 .*

*bSheep density: number of sheep/km2*

*c Pig density: number of pigs/km2 .*

*dHuman density: number of humans/km2 .*

each epidemic. The median number of IPs on days 7, 14, and 21 for ISP was 6, 9, and 11 (respectively) compared with 3, 5, and 5 for AADIS (**Tables 2** and **3**). Similarly, the number of traces generated by ISP in the early phase of each epidemic was higher than those generated by AADIS. The median number of traces generated by days 7, 14, and 21 for ISP was 8, 11, and 12 (respectively) compared with 2, 3, and 4 for AADIS. There are three possible explanations for these findings: (1) differences in characteristics of the countries and/or study regions and incursion scenarios used for each model; (2) differences in model parameterization, resulting in different probabilities of farm-to-farm transmission of virus; and (3) differences in model design (in ISP the probabilities of transmission vary according to farm type but not farm size whereas in AADIS both farm size and farm type influence probabilities of transmission).

Outbreak durations for the two models were similar: a median of 43 (minimum 16, maximum 365) days for AADIS compared with a median of 43 (minimum 21, maximum 263) for ISP. The size of the AUC was substantially lower for the AADIS simulations. The median AUC for the AADIS simulations was 680 km2 (minimum 300, maximum 29,953) compared with 1176 km2 (minimum 316, maximum 12,815) for ISP.

#### Linear Regression

Regression coefficients and their standard errors for the linear regression models of the total number of infected places, outbreak duration, and the total AUC for the AADIS and ISP models of FMD are provided in **Table 5**. **Table 6** provides details of the goodness of fit (*R*2 ) for each of the linear regression models developed for Australia and New Zealand. A consistent pattern was observed with the goodness of fit of the models improving from days 7 to 14 to 21 for all outcome variables for both the Australian and New Zealand data sets.

Positive and negative predictive values for "large" or "small" outbreaks (for the total number of IPs and total AUC) or "short" or "long" outbreaks (for outbreak duration) for AADIS and ISP

*EDR, estimated dissemination ratio.*

TABLE 6 | Goodness-of-fit statistics (*R*2 ) for each of the linear regression models for the total number of infected places, outbreak duration, and area under control using days 7, 14, and 21 explanatory variables for the AADIS and InterSpread Plus models of FMD.

are shown in **Table 7**. The proportions of correctly classified outbreaks ranged from 0.88 to 0.97 for AADIS and 0.79 to 0.92 for ISP.

#### Regression Trees

Classification and regression tree analyses were carried out to identify factors associated with the total number of IPs, outbreak duration, and total AUC. Similar to the approach used for the linear regression analyses, three sets of explanatory variables were used: those at day 7 post-detection, day 14 post-detection, and day 21 post-detection. Using these three sets of explanatory variables with each of the three outcome variables and both the Australian and New Zealand data sets resulted in 18 CART analyses in total. BRT models using the same explanatory variables and the same outcome variables were developed using the Australian and New Zealand data.

The CART for the predicted total number of IPs using day 14 explanatory variables for the Australian and New Zealand data are shown in **Figures 2** and **3**, respectively. For both the AADIS and ISP models, the number of IPs at day 14 had the greatest influence on the total number IPs. For the AADIS model, in TABLE 7 | Positive and negative predictive values and the proportion of outbreaks correctly classified as large or small (or short or long) using the day 14 linear regression model for the AADIS and InterSpread Plus simulated FMD outbreaks.

*a Used to classify outbreak as small or large.* addition to the number of IPs identified at day 14, the total AUC at day 14 and cattle density influenced the total number of IPs.

The five most influential explanatory variables (and their weights) from the BRT models for the total number of IPs, outbreak duration, and the total AUC for the Australian and New Zealand data are listed in **Table 8**. Consistent with the CART analyses, for both AADIS and ISP, the number of IPs identified at day 14 was associated with each of the three outcome variables. While for ISP the number of IPs identified at day 14 had the highest weight for each of the three outcomes, the total number of outbreak clusters identified at day 14 had the greatest weight as a predictor of the total AUC for AADIS. For AADIS, the density of cattle was associated with each of the three outcome variables, albeit with a relatively low regression weight in each model (10.2, 15.4, and 0.1 for the total number of IPs, outbreak duration, and the total AUC, respectively).

The predictive ability of each of the day 14 boosted regression models was assessed by calculating the positive and negative predictive values for each model (**Table 9**), similar to the

TABLE 8 | Identified explanatory variables (*n* **=** 5) and their weights (in brackets) for the boosted regression tree model of first 14 day predictors of area under control, the total number of infected places, and outbreak duration for the AADIS and InterSpread Plus models of FMD.

approach taken for the linear regression models. Overall, the BRT models were able to correctly classify simulated outbreaks as either large or small (for the total number of IPs and total AUC) or short or long (for outbreak duration) with the proportion of correctly classified outbreaks ranging from 0.82 to 0.96 for AADIS and 0.77 to 0.93 for ISP. In general, negative predictive values for the BRT models were greater than the positive predictive values.

#### DISCUSSION

During a disease outbreak, decisions on control are often made under significant uncertainty and in conditions that are continually evolving. Resources are often limited and will influence the effectiveness of disease control efforts. Experience overseas suggests that resource and logistical issues are critical considerations when evaluating disease control strategies (39–41). Vaccination is increasingly being recognized as an important tool to assist in containing and eradicating FMD outbreaks (6–8, 42, 43). Garner et al. Early Decision Indicators for FMD Control

TABLE 9 | Positive and negative predictive values and the proportion of outbreaks correctly classified as large or small (or short or long) using the day 14 boosted regression tree model for the AADIS and InterSpread Plus simulated FMD outbreaks.

*a Used to classify outbreak as small or large.*

Vaccination has been shown to be most effective in situations where disease is spreading rapidly or resources are inadequate to maintain effective stamping out (4). A number of studies have shown that vaccination is more beneficial when used early in an outbreak (8, 44, 45).

Although vaccination can be an important tool to control FMD, it will make achieving recognition of FMD-free status more difficult – keeping vaccinated animals in the population will delay the period until FMD-free status is regained under the World Organization for Animal Health (OIE) guidelines and add additional complications to the post-outbreak surveillance program (46). Shifting attitudes to vaccination among the international veterinary community means that it is no longer viewed as a measure of last resort. In Australia and New Zealand, vaccination will be given consideration as a potential additional measure (alongside stamping out) from day one of any FMD eradication response. However, given the complications and costs associated with implementing a vaccination strategy, it would only be used if authorities consider that it would be beneficial in managing the outbreak (19, 26, 47). A decision to vaccinate early in the outbreak may result in situations where it was not actually required and have consequent implications for post-outbreak surveillance, management of vaccinated animals, and regaining FMD-free status and access to markets. Conversely, not using vaccination in some situations may lead to larger and longer outbreaks, increased control costs, and greater on-going impacts on industry and local communities.

Although Australia and New Zealand have developed frameworks to support decision-making on FMD control (19, 48), these are qualitative and subjective. We reason that it would be useful if disease managers could identify early in an outbreak those situations that are likely to progress to "large" outbreaks and for which additional measures like vaccination are likely to be beneficial. In this context, measurable parameters such as the number of IPs, numbers of traced premises and/or farms under surveillance, and estimated rates of spread might be useful indicators of the potential severity of an outbreak.

The overarching aim of this project was to identify factors that could be used to predict the total number of IPs, outbreak duration, and the total AUC. Here, "factors" refers to characteristics of the physical environment in which an FMD incursion first occurs (e.g., farm density, animal density, human population density) or characteristics of the outbreak itself (e.g., the number of IPs reported at a given point in time post first detection). We were particularly interested in how robust the findings were to outbreaks in different settings. For this study, we used a wide range of FMD incursions in terms of location, production systems and seed farm type, and time to first detection (determined probabilistically). These outbreaks were simulated in two countries using two independent modeling platforms.

It is reassuring for animal health authorities that, in both countries, the simulated FMD outbreaks tended to be small and readily able to be contained and eradicated with available resources. For both countries, median outbreak durations were around 6 weeks. This finding assumes that FMD is reported relatively quickly and resources are adequate to implement effective control programs. For Australia, the median time from first introduction to reporting was 17 days (range 9–89), and for New Zealand, the median time to detection was 13 days. A previous Australian study found considerable regional variability in the probability that an individual infected farm would report suspect FMD (24, 49). Recent experience of outbreaks of FMD in non-endemic countries indicate that it can take up to 3 weeks after introduction of the virus to the primary farm before the disease is recognized (40, 50–52). However, early detection does not necessarily mean that an outbreak will be small. A total of 3.4% of the 10,000 outbreaks of FMD in Australia that were simulated in this study had more than 100 IPs and 7.2% of the 10,000 outbreaks lasted longer than 90 days. For New Zealand, there was a 7.2% probability of an outbreak involving more than 100 IPs and an 8.6% probability of an outbreak lasting more than 90 days.

The key objective of this study was to test whether information known or available to disease managers early in an FMD outbreak could be used to predict the severity of the epidemic outcome. Epidemic outcome was defined in terms of the total number of IPs, outbreak duration, and the total geographic AUC. While FFO and FFS have been shown to correlate with epidemic size (16), it was recognized that it would be more useful to consider a broader range of times than just 14 days. Accordingly, three time points were considered: 7, 14, and 21 days into the control program. A range of potential explanatory variables were tested using different analytical approaches, including linear regression, CART, and BRT analyses.

Although there was some variability between the different analyses and between countries, the cumulative number of IPs at specified time points early in the outbreak were consistently found to be strongly associated with the final number of IPs and the duration of an outbreak. It was possible to build relatively simple linear regression models for predicting the magnitude and duration of simulated FMD outbreaks that fitted both the Australian and New Zealand data (see **Table 5**). *R*2 values as a measure of goodness of fit ranged from 0.3 to 0.9 depending on time point, outbreak variable, and country (**Table 6**). A consistent pattern was observed, with the fit of the models improving from days 7 to 14 to 21 for all dependent variables and for both data sets (Australia and New Zealand). The total AUC had the highest predictability and duration of an outbreak the lowest. In this study, we found that the number of IPs occurring up to a given time point provided the most predictive power for both size (total IPs) and outbreak duration. This confirms previous findings by Hutber et al. (15) and Halasa et al. (16). The AUC at a given time point was most predictive of the total AUC.

These findings were confirmed in the CART and BRT analyses. Consistency between the different approaches helps build confidence that the criteria identified are relevant to response decision-making. CART techniques are a useful alternative as they provide a visual decision tree output that is intuitive and likely to be well received by those not familiar with statistical analysis (see **Figures 2** and **3**). The tree diagrams produced in a CART analysis are consistent with clinical reasoning used by animal health professionals and can help to structure explanations of prediction. Compared with regression-based approaches, an advantage of a CART analysis is that it can accommodate non-linear relationships between an outcome variable and a set of explanatory variables as well as missing data.

Boosted regression trees have the advantages of being able to handle a range of explanatory variable types, not requiring any data transformations, and being able to account for complex, non-linear relationships (35). BRTs are better-able to describe linear relationships and are more robust in terms of predictive accuracy, although interpretability suffers as a result. CARTs and BRTs are complementary. CARTs are relatively simple and provide readily interpretable output; BRTs are more complex and robust, but with reduced interpretability. The BRTs for both countries had good predictive ability when the total number of IPs was less than 100. When the total number of IPs was greater than 100, the BRT analyses tended to under predict total IP numbers.

Although it is informative to build statistical models to summarize factors influencing outputs from complex simulation models of FMD, for disease managers, the key issue is how this information can be used to support decision-making. From a disease manager's perspective, it is useful to consider how good the models are at predicting small and large outbreaks. To do this, it is necessary to make some judgment calls about what constitutes a "large outbreak." It is difficult in advance to reach agreement on what are acceptable benchmarks in terms of eradicating FMD, as this will be influenced by the time and location of an outbreak, availability of resources, etc. Accordingly, we looked at a series of arbitrary "cut points" for classifying outcomes into small and large (or long and short) outbreaks. Model sensitivity, specificity, and positive and negative predictive values were calculated using these cut-points. In general, the linear regression models were very good at predicting when an outbreak would be small or short; the positive predictive values varied from 0.85 to 0.98 meaning that a small outbreak was correctly predicted between 85 and 98% of the time. It should also be noted that having predicted a small outbreak at day 14 (which would probably mean that a decision to vaccinate would not be made), this decision could be revisited at a later time in the outbreak when more information was available. Incorrectly predicting a large outbreak and using vaccine when it is not actually required will have trade implications and increase outbreak costs. The models were less accurate at predicting a large or long outbreak with the negative predictive values for outbreak duration exceeding 90 days being as low as 0.52 for the models of FMD in New Zealand. The negative predictive values for the total number of IPs and the total AUC were better ranging from 0.77 to 0.91 for both AADIS and ISP.

The BRT models were able to correctly classify simulated outbreaks as either large or small with the proportion of correctly classified outbreaks ranging from 0.77 to 0.96. Negative predictive values tended to be higher than the positive predictive values for the BRT models.

In conclusion, this study shows that based on simulated FMD outbreak data relatively simple metrics available at 1–3 weeks into the control program can be used to predict the size of an FMD outbreak under Australian and New Zealand conditions and provide a basis for making decisions on the use of vaccination as a control measure. It should be noted that the simulation modeling analyses carried out for this study focused on introduction of FMD into the areas considered to be at higher risk of disease entry and dissemination in Australia and New Zealand (23). The results need further validation with modeling data generated from other areas of these countries. Finally, it should be recognized that in the absence of FMD outbreaks in Australia and New Zealand, this study has fitted statistical models to simulated, not real outbreak data. Although the modeling teams have been careful to parameterize the respective models as realistically as possible, it is inevitable that assumptions and extrapolations from overseas experience have had to have been made. These considerations need to be taken into account when using the findings from this study.

#### AUTHOR CONTRIBUTIONS

All the authors made substantial contributions to the work, through concept and design (MG, TR, RS, SR, TK, and PH); data generation and modeling (MG, TR, RS, and RB); statistical analyses (IE, MS, and RS); or interpretation (MG, IE, MS, TR, and TK). MG, IE, and MS had prime responsibility for drafting and revision. MG, TR, and TK have approved the final version for publication and are in agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

#### FUNDING

This work was funded by the Centre of Excellence for Biosecurity Risk Analysis (CEBRA) Project 1404D and supported by the Australian Department of Agriculture and Water Resources and the New Zealand Ministry for Primary Industries.

#### REFERENCES

disease in animal populations. *Prev Vet Med* (2013) 109:10–24. doi:10.1016/j. prevetmed.2012.08.015

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer PW and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

*Copyright © 2016 Garner, East, Stevenson, Sanson, Rawdon, Bradhurst, Roche, Van Ha and Kompas. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Semiquantitative Decision Tools for FMD Emergency Vaccination Informed by Field Observations and Simulated Outbreak Data

*Preben William Willeberg1 \*, Mohammad AlKhamis2,3, Anette Boklund1 , Andres M. Perez3 , Claes Enøe1 and Tariq Halasa1*

*1Department of Diagnostic and Scientific Advice, National Veterinary Institute, Technical University of Denmark, Copenhagen, Denmark, 2Environment and Life Sciences Research Center, Kuwait Institute for Scientific Research, Kuwait City, Kuwait, 3Department of Veterinary Population Medicine, College of Veterinary Medicine, University of Minnesota, St. Paul, USA*

#### *Edited by:*

*Alejandra Victoria Capozzo, Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Argentina*

#### *Reviewed by:*

*Fedor Korennoy, Federal Center for Animal Health (FGBI ARRIAH), Russia Francois Frederick Maree, Agricultural Research Council, South Africa*

*\*Correspondence:*

*Preben William Willeberg<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Veterinary Epidemiology and Economics, a section of the journal Frontiers in Veterinary Science*

*Received: 30 November 2016 Accepted: 09 March 2017 Published: 27 March 2017*

#### *Citation:*

*Willeberg PW, AlKhamis M, Boklund A, Perez AM, Enøe C and Halasa T (2017) Semiquantitative Decision Tools for FMD Emergency Vaccination Informed by Field Observations and Simulated Outbreak Data. Front. Vet. Sci. 4:43. doi: 10.3389/fvets.2017.00043*