Drivers for Livestock-Associated Methicillin-Resistant Staphylococcus Aureus Spread Among Danish Pig Herds - A Simulation Study

To gain insight into the rapid increase in the number of livestock-associated Methicillin-resistant Staphylococcus aureus (LA-MRSA)-positive herds in Denmark, we developed an individual-based Monte Carlo simulation model. We aimed to assess whether transmission of LA-MRSA via pig movements could explain the observed increase in the number of positive herds in Denmark, and to evaluate the effect of other between-herd transmission mechanisms. Pig movements alone were not sufficient to mimic the observed increase in LA-MRSA-positive herds in Denmark in any of the modelled scenarios. The model identified three factors that played important roles in the between-herd spread of LA-MRSA: (1) the within-herd dynamics, (2) the frequency and effectiveness of indirect transmissions, and (3) unexplainable introduction of LA-MRSA to swine herds. These factors can act as starting points for the development of LA-MRSA control programs in pig herds in order to limit the risk of its transmission to humans.


Data preparation
Herd information. Holding information data included the holding and herd identification numbers.
Several herds could be owned by the same farmer and therefore share an individual holding number, yet have different herd identification numbers. To keep the herds anonymous, new herd identification numbers were generated, and the herd was used as the unit of interest. Herds owned by the same farmer were marked to track the connection between them. The registered numbers of sows, weaners and finishers per year, as well as the herd type were extracted from the CHR.
The VetStat database records herd-level information on prescription-only drugs 1 . These data were used to add an antimicrobial usage index per year. The index was set to 1 if at least one prescription of tetracycline and/or beta-lactam was reported in the VetStat register.
In addition, the number of registered movements to abattoirs for each herd per year was calculated based on the movement data and this value divided by 365, giving a herd-specific that was used as a parameter to model indirect contacts related to abattoir movements.
Our datasets lacked information on UTM coordinates for some of the herds. This information was updated using a website to calculate the coordinates based on the address of the holdings provided on the CHR website 2 (www.geoplaner.com/).
All herds were categorised according to: (1) herd type (breeding and multiplier herds, production herds, weaner herds, organic and free-range pig herds and hobby herds), (2) the proportion of sows, weaners and finishers registered in the CHR 2 , and (3) the type of production (Table S1-S2).
No sows, weaners or finishers were registered in 181 herds in 2007 and 844 herds in 2014 3 , but these herds had registered in-coming or out-going movements, and could therefore be considered active. The number of animals was estimated for these herds. For each herd type, the distribution of herds (with registered sows, weaners or finishers) was determined in each combination of categories 1-7 and A-C (Tables S1-S2). In addition, the average number of sows, weaners and finishers registered in these herds in the CHR was calculated for each category combination. To estimate the missing number of animals, a category combination was randomly assigned to herds that had not registered any animals.
Based on this combination, the number of sows, weaners and/or finishers was calculated using an exponential distribution with lambda given by the average number of sows, weaners and finishers.
Movement data. Information on the movements of swine was available from the movement database, including the holding and herd identification numbers for both sending and receiving herds, the date of the pig movements and the number of pigs moved. No information was available on the age group (sows, weaners or finishers) of pigs moved.
The age group of pigs moved out of the sending herd and into the receiving herd was estimated based on the herd categories 1-7 and A-C (Tables S1-S4). If no weaners were registered, but were thought to have been moved out of a herd, they were assumed to have been moved directly from the sow section.

Modelling disease spread within a herd
Environment-related recurrence. The probability decreases exponentially over time as follows: where is defined as the time difference between the current simulation day and the day when LA-MRSA died out in the respective herd.
Within-herd dynamics. For each individual susceptible animal, the probability of infection at each time step resulting from positive animals within compartment was assumed to be densitydependent and was calculated based on the number of LA-MRSA-positive pigs and the total number of pigs within the individual compartment as follows: where is the daily transmission rate of the infection for each compartment , and The probability of infection resulting from positive animals from another compartment was calculated as: where is the between-compartment transmission rate. The total probability for each individual in compartment was therefore: In order to slow the within-herd spread, a time-dependent scaling function was introduced: where is the simulation time since the introduction of LA-MRSA in the herd, and and are readin parameters defining the steepness and the midpoint of the scaling function. The within-compartment transmission rates were divided by the scaling function. The time-shifted transmission rates led to a similar but delayed level of within-herd prevalence (Figures S1-S5).
Mimicking the within-herd dynamics in 100 large production herds for 365 days using the transmission rates adapted by Broens et al. (2012) 4 led to a median within-herd prevalence of 65% for herds using high-risk antibiotics and 59% for herds that did not use high-risk antibiotics (Figures S1-S4).
Spread via indirect contacts. The transmission between two herds and with the same owner was modelled using equation (2) Figure S1. Development of the within-herd prevalence in large production herds in the first year after LA-MRSA initialisation in 100 production herds (1) using high-risk antibiotics, and (2) not using high-risk antibiotics for (a) transmission rates adapted from Broens et al. 4 (Table S5), and (b) time-shifted transmission rates. The black dotted line represents the median of 500 iterations and the grey area spans the 95% confidence interval. Figure S2. Development of the within-herd prevalence in the sow compartment in the first year after LA-MRSA initialisation in 100 production herds (1) using high-risk antibiotics, and (2) not using high-risk antibiotics for (a) transmission rates adapted from Broens et al. 4 (Table S5), and (b) time-shifted transmission rates. The black line represents the median of 500 iterations and the grey area spans the 95% confidence interval.         Table S3. Assumed age group (s = sows, w = weaners, f = finishers) of pigs moved out of the sending herd / in to the receiving herd. The letter x is used when the movement type is dependent on a defined threshold (Table  S4). Below or equal to this threshold, finishers are assumed to be moved out of the sending herd and sows are assumed to be moved in to the receiving herd. Above this threshold, it is assumed that weaners/finishers were moved in to / out of the sending / receiving herd. 1A

TABLES
w/f w/f w/f w/f w/f w/f w/s w/f w/s w/f w/w Table S4. Assumed thresholds for setting the age group of movements if the movement type is given as x in Table S3. The thresholds are dependent on the registered herd type and for production, organic and hobby herds on the size of the herd. In small herds, homogeneous mixing of pigs is assumed, whereas separate compartments for the three age groups (sows, weaners, and finishers) were modelled for large herds. The number of pigs assumed to constitute a small herd was set to 200 animals. Table S5. Assumed values for a PERT distribution to define herd-specific cure rates and transmission rates based on the use of high-risk antibiotics, adapted by Broens et al. 4 The most likely (mode) and assumed minimum and maximum values for the PERT distributions were calculated based on values for R 0 (and their 95% CI) resulting from multivariable analysis of Dutch data, with 10.3 days taken to be the duration of the infectious period 5 . PERT distributions were defined as transformation of the Beta distribution with minimum (min), maximum (max) and most likely value (mode) and a mean .