Modelling spread and surveillance of Mycobacterium avium subsp. paratuberculosis in the Swedish cattle trade network.

To monitor a state of disease freedom and to ensure a timely detection of new introductions of disease, surveillance programmes need be evaluated prior to implementation. We present a strategy to evaluate surveillance of Mycobacterium avium subsp. paratuberculosis (MAP) using simulated testing of bulk milk in an infectious disease spread model. MAP is a globally distributed, chronic infectious disease with substantial animal health impact. Designing surveillance for this disease poses specific challenges because methods for surveillance evaluation have focused on estimating surveillance system sensitivity and probability of freedom from disease and do not account for spread of disease or complex and changing population structure over long periods. The aims of the study were to 1. define a model that describes the spread of MAP within and between Swedish herds; 2. define a method for simulation of imperfect diagnostic testing in this framework; 3. to compare surveillance strategies to support surveillance design choices. The results illustrate how this approach can be used to identify differences between the probability of detecting disease in the population based on choices of the number of herds sampled and the use of risk-based or random selection of these herds. The approach was also used to assess surveillance to detect introduction of disease and to detect a very low prevalence endemic state. The use of bulk milk sampling was determined to be an effective method to detect MAP in the population with as few as 500 herds tested per year if the herd-level prevalence was 0.2 %. However, detection of point introductions in the population was unlikely in the 13-year simulation period even if as many as 2000 herds were tested per year. Interestingly, the use of a risk-based selection strategy was found to be a disadvantage to detect MAP given the modelled disease dynamics.

Introduction

Paratuberculosis is a chronic infectious disease of ruminants caused by Mycobacterium avium subsp. paratuberculosis (MAP). It is globally distributed and has a substantial animal health and economic impact in dairy and beef farming (Whittington et al., 2019) and potentially a threat to human health (Waddell et al., 2015). Sweden has a unique low prevalence and possibly even an absence of MAP in cattle herds (Frössling et al., 2013; Lindberg, 2019; Whittington et al., 2019). The strategy to control, eradicate and prevent the introduction of MAP in Sweden was established early in the 20th century (Engvall et al., 1994). Furthermore, Swedish legislation enacted in 1952 made the notification and investigation of all clinical suspicions of paratuberculosis mandatory in all species. Trace-back investigations and eradication measures, including whole-herd depopulation have subsequently been completed upon detection (Engvall et al., 1994; Sternberg and Viske, 2003; Frössling et al., 2013). Historically, these detected cases have been found to be linked to the import of animals (Lindberg, 2019). Sweden maintains a database of all cattle, each unique holding identifier in the database corresponds to a location where cattle are kept, for example, a farm building or a pasture. In the same way as in Widgren et al. (2016), these are hereafter jointly defined as holdings. One of the main components of the surveillance system for MAP in Swedish cattle is a specific voluntary surveillance programme that covers mainly pedigree beef holdings and is run by the Swedish Animal Health Service (Frössling et al., 2013; Lindberg, 2019). This voluntary programme started in 1998 and affiliated holdings have agreed to only purchase animals from holdings with the same certified status.

Due to imperfect testing, it is not directly possible to prove the absence of a disease in a population. It is only possible to apply tests to the population to detect a disease, and over time, consistently negative test results contribute to evidence in favour of freedom. However, the collection of such evidence of MAP freedom is resource intensive because of the slow clinical progress of the disease, slow spread and low diagnostic sensitivity in latently infected animals. The contribution of a surveillance effort to the body of evidence supporting freedom from disease can be estimated using the sampling frequency, distribution of disease, diagnostic performance and a defined lower detection threshold (Martin et al., 2007). Frössling et al. (2013) used stochastic scenario-tree modelling to estimate the probability of freedom from MAP in Swedish cattle. This method can be used to combine information from several sources, to incorporate differences in risk and probability of detection for different population strata, and also to account for the probability of disease introduction. This approach is focused on estimating surveillance system sensitivity and probability of freedom from disease. However, when designing a surveillance programme, other performance parameters may also be relevant such as the time to detection of an introduction of disease or the number of affected herds at the time of detection. Additionally, scenario-tree modelling does not capture the contacts between holdings that change over time, or the spread between risk groups which are described in systematically recorded animal movements between holdings. For example, during the period 2005–2017, the total number of holdings with cattle in Sweden decreased from 26,179 to 16,674 (Statistics Sweden, 2008, 2018). Approximately 60 % of the adult cows were dairy cows, and 40 % were suckler cows, but due to differences in holding size, suckler holdings (cow-calf operation) were more common than dairy holdings. In addition, there were a small number of specialised beef holdings that rear calves and youngstock for slaughter and these calves were mainly purchased from dairy holdings. Such dramatic changes in population structure can affect the contact patterns over time and hence disease dynamics, which in turn has an impact on surveillance performance.

In a recent study, VanderWaal et al. (2017) presented a strategy to evaluate surveillance efficacy by simulating disease spread, surveillance and control in a dynamic temporal network model. This strategy could be applied in an infectious disease modelling framework that uses recorded data (births, deaths, movement of animals between holdings) to drive population demographic and disease spread between herds over time. To formalise the simulation of disease spread and surveillance on a temporal network would allow for the optimisation of surveillance to fit the diverse disease spread scenarios and surveillance targets. Also, these targets can be expanded from surveillance system sensitivity or probability of freedom to any parameter that is measurable from a disease spread model such as the time to detection and the extent of an outbreak at the time of detection.

The first aim of this study was to define a model that describes the spread of MAP within and between herds to investigate potential infection patterns and network features relevant for disease spread, given the existing Swedish cattle population and different disease introduction scenarios. The second aim was to define a method for simulation of imperfect diagnostic testing in a disease spread model. The third aim was to compare surveillance strategies and support surveillance design choices and conclusions about MAP occurrence or absence in Sweden.

Material and methods

In this study, we used the modelling framework SimInf, described by Widgren et al. (2019) to simultaneously model within-holding and between-holding spread of MAP in the Swedish cattle population over the study period: 2005–2017. To include the specific aspects of within-holding spread of MAP, a within-holding model of MAP was adapted from a previous model described by Smith et al. (2015). Animal events, including movements, birth, and ageing, were included in the model by using real recorded animal life-event data from the Swedish national cattle database in the simulations. Different surveillance alternatives based on testing of bulk milk were then applied in the model.

Input data

To incorporate the complete Swedish cattle population and livestock movements in the modelling, all records of birth, death or movements of cattle in the period 2005–2017 were retrieved from the Swedish national cattle database. The data included 40,349 unique holding identifiers and 15,237,843 recorded births, deaths or movements which are reported digitally at the individual animal level within 7 days of the actual event to the Swedish Board of Agriculture. The change in number of holdings, animals, births, deaths and movements over time and season has been previously described in detail for the Swedish cattle industry (Widgren et al., 2016). The 10th, 25th, 50th, 75th and 90th percentiles of herd size in June 2009 were: 3, 9, 28, 76 and 151 respectively. This register data was cleaned in that same way as previously described (Widgren et al., 2016). To enable analyses stratified by production type, a list of dairy holdings kept by the Swedish dairy association (Växa Sverige) was used (4715 holdings). Any holding in that list was defined as dairy and others were considered to be non-dairy. To describe the age distribution of animals in the dairy sector, in the average dairy herd on 2015-06-15: 39 % of female animals were < 2 years of age, 15 % were between 2 and 3 years, 12 % between 3 and 4 years, 9% between 4 and 5 years and 24 % were ≥ 5 years of age. A list of members sampled in the voluntary surveillance programme for MAP in 2008 (221 holdings) were obtained from the Swedish Animal Health Service to study the spread of MAP between members of the programme and those holdings outside the programme. Of these holdings, 218 were non-dairy and 3 were dairy. There is no control programme, such as, test and cull for MAP in Sweden and MAP vaccine is not used in the population. For investigation of potential risk-based surveillance strategies, the number of unique holdings that moved animals directly into each holding was determined for all holdings over the period 2005–2017. This metric is known as indegree (Wasserman and Faust, 1994).

Within-holding model

The within-holding spread of MAP, was modelled by adapting a model previously presented by Smith et al. (2015). This is an age-stratified (calves, heifers and cows) compartment model of MAP consisting of susceptible and transiently or latently infected high or low shedders (Fig. 1 ). The transitions between the compartments were modelled as a continuous-time discrete state Markov chain using the Gillespie’s direct method (Gillespie, 1977). A more detailed description of the current implementation including all propensity functions and their justification is presented in the Appendix 1 - Supplementary material. The within-holding model was implemented using the SimInf package version 6.1.0 in R version 3.5.1 (R Core Team, 2018).

Fig. 1

A schematic representation of the model of MAP. For simplicity, births, deaths and trading of animals between holdings (external transfers) are not included in the diagram. All births go to S∅ and deaths and external transfers may happen from all nodes in the model. Calves are defined as animals less than 1 year of age, heifers between 1 and 2 years of age and cows older than 2 years. The propensities of all stochastic processes in the model are described in detail in Appendix 1 - Supplementary material. The transitions from S∅ to S1, T1H, or T1L in calves was proportional to infection in cows to simulate vertical transmission. Heifers could be susceptible (S2), transient (T2H, or T2L) or latently infected (E2H or E1L). Cows could additionally enter the shedding state either high (L2H, H3) or low (L2L). H and L in the compartment names indicates hight and low shedder respectively. Note that cows newly infected as adults could only join the low shedding group. Additional compartments were included in the model to monitor the test results when simulating surveillance (T∅, Tpos, Tneg).

Initialisation of the within-holding model was done by moving animals in each age category from the susceptible compartment to one of the infected states. The number of animals to be moved to an infected state was calculated by rounding 2% of the number of animals in each age category up to the nearest integer. The initial infection was chosen to be small but in multiple age-categories of the holding to reduce the chance of stochastic extinction of the infection. To explore the outcome from using this initialisation scheme, 10,000 trajectories (one trajectory is the outcome from one simulation) were generated for each of: a small (12 animals), a medium (32 animals) and a large (76 animals) holding, respectively. These holding sizes corresponded to the 25th, 50th and 75th percentile of the holding sizes in Sweden in 2005. The median and 95 % credibility interval around the number of infected animals in each holding over time was plotted, where the 95 % credibility interval was defined as the range between the 2.5 % and 97.5 % percentiles at each time step and calculated by the default quantile method in R version 3.5.1.

National-scale model

In the national-scale model, the within-holding spread model was run within the network of Swedish cattle holdings to model between-holding spread of MAP. In this implementation of the model within a network of holdings, the births, ageing, death and trade of animals were controlled by the real recorded animal births, deaths and movement events and has been previously described (Widgren et al., 2019). Events for each day and holding were applied following the stochastic part of the model which moved animals between the states of infection within a holding. The only means of spreading infection between holdings in this model was via animal movements. The national scale model was run starting with initialisation in a single holding, completing 12 trajectories for each of the 28,502 holdings with animals in 2005.

The number of secondary holding infections were recorded to determine the impact of seeding the infection in different holding types on the disease spread. The model results were summarised to determine the spread between dairy and non-dairy holdings and within dairy and non-dairy respectively. To assess how the voluntary surveillance programme member and non-member holdings could contribute to spread dynamics of MAP respectively, the results were summarised to determine spread between programme-member and non-member holdings. The median, mean and 95 % credibility interval around the number of infected holdings was plotted.

Surveillance

Simulation of disease testing was implemented in the model to enable the investigation of the performance of varying surveillance efforts. The SimInf framework allows for the inclusion of pre-scheduled events to change the counts in compartments at specific time-points. Three compartments were added to the model to perform disease testing: 1) an indicator compartment T ∅, to switch-on the sampling in a holding at a specific point in time 2) T to indicate a positive test result, and 3) T, to indicate a negative test result (Fig. 1). The indicator compartment T ∅ was switched-on by injecting a sampling event into the model to change T ∅ from 0 to 1. This triggered two possible paths, either T ∅ was reclassified to T with a rate proportional to the diagnostic test sensitivity of the bulk milk ELISA (SeBM) if any infected adult animals were present in the holding, or to T, proportional to one minus the test sensitivity. The diagnostic test specificity was assumed to be 1.

The diagnostic test sensitivity of the bulk milk ELISA (SeBM) can be interpreted as the probability of the bulk milk ELISA test being positive given a single infected adult animal in the holding. SEBM was assumed to be dependent on the prevalence of infected cows in the holding since this would affect the concentration of antibodies in the bulk milk tank. To estimate this relationship, 400 data points were extracted from Fig. 3 of Nielsen and Toft (2014), who reported the result of bulk milk testing from holdings with varying apparent within-holding prevalence. The cut-off value for the sample to positive (S/P) ratio of the test to declare a sample positive was set to 10 %. This cut-off value was chosen rather than the 15 % used in previous work in order to improve the test sensitivity. A logistic regression was used to estimate the association between the bulk milk ELISA binary test result (outcome) and the log2 transform of the apparent within-holding prevalence (explanatory variable). Then this model was used in the simulations to calculate SeBM from the simulated apparent within-holding prevalence each time a test was to be applied. The apparent prevalence was calculated by multiplying the true prevalence in the simulator with the diagnostic test sensitivity for an animal test (SeI) estimated by Nielsen et al. (2013) and adjusted to the age distribution in a Swedish herd. This resulted in an SeI of 62.7 % being used in the simulations. The methods employed to estimate SeI and the implementation of surveillance in the model are described in more detail in Appendix 1 - Supplementary material.

Fig. 3

The median (solid black line), average (dashed black line) and 95 % credibility interval (grey area) of the number of secondary infected non-dairy or dairy holdings that remain infected over time after an introduction of MAP infection into a single either non-dairy or dairy holding. Simulations included all Swedish cattle farms 2005–2017. Note, the median number of secondary cases is 0 in all 4 graphs.

Surveillance for MAP using bulk milk testing was evaluated by applying testing to 500, 1000, 1500, or 2000 holdings per year and targeting holdings for testing by random or risk-based selection. Risk-based selection was either weighted by normalised 13-year log2(indegree) or weighted by normalised 13-year indegree. To capture different aspects of surveillance performance, the following outcomes were measured for each surveillance strategy: 1) The surveillance system sensitivity i.e., the proportion of detected outbreaks; 2) the time to outbreak detection and; 3) the size of the outbreak at the time of detection. Survival analysis and the Cox proportional hazards model were used to compare the cumulative hazard of detection (the expected number of detections) between surveillance strategies over time (Dardis, 2018). This regression model included a categorical variable for number of sampled holdings, with 500 as the reference group and the holding-selection method with random selection as the reference group. Trajectories were censored at the first occurrence of either: outbreak die out, outbreak detection or the end of the 13-year simulation period. Hazard ratios significantly different from 1 at a confidence level of 95 % were considered to indicate two surveillance strategies with significantly different hazard of detection. The Cox proportional hazards regression approach was used to interpret both time to detection and the surveillance system sensitivity together as these parameters of a metrics of surveillance are related to one another. The assumption of proportional hazards of the Cox model was evaluated by visualising the Schoenfeld residuals over predictor and time. The model fit was evaluated by plotting the Cox-Snell residuals against the cumulative hazard.

To evaluate the ability of the surveillance system to detect point introductions while disease spread was limited to a small part of the population, trajectories were generated from the national-scale model with infection seeded at a single holding. For each holding with animals in 2005, 12 trajectories of the model were run seeding the infection in 1 herd, i.e. each of the 12 surveillance strategies was applied once to each starting herd scenario, producing a total of 342,024 trajectories. Secondly, to evaluate the efficacy of the proposed surveillance strategies to detect a low prevalence, the national-scale model was initialised with infection seeded in a random sample of 0.2 % of holdings (57 holdings in 2005). This seeding was repeated 1000 times for each surveillance strategy, producing a total of 12,000 trajectories. This between-holding prevalence was chosen to reflect a lower detection limit of 0.2 % between-holding prevalence, previously used as a surveillance design target by Frössling et al. (2013).

Results

Spread within a single closed holding

Results from the three scenarios used to investigate the behaviour of the model in a single holding and assess the efficacy of seeding are presented in Fig. 2 . For these very small holding sizes and limited within-holding prevalence of MAP, the model tended toward extinction of the disease. The extinction times were similar in the various holding sizes tested, and were consistently longer than the 13-year study period used in the national scale model.

Fig. 2

Progression of the number of animals infected in a holding after seeding infection in a holding with a total holding size of 12 (6 cows, 3 heifers, 3 calves), 32 (16 cows, 8 heifer, 8 calves), and 76 (38 cows, 19 heifers, 19 calves) animals (left to right), representing a small (25th percentile) median and large (75th percentile) holding in Sweden. The solid black line indicates the median number of infected animals over time and the grey areas the 95 % credibility interval. Seeding of infection in the holdings was completed in the same way as other models in the current study by rounding up 2% of the animals stratified by age group to be moved to one of the infected states.

Spread between subgroups in the Swedish cattle population

Based on the model and study period, the median number of secondary infections after a single introduction in any type of cattle holding was <1. Stratification of these results by introduction in the dairy and non-dairy sectors (Fig. 3) showed that spread was relatively more common from dairy holdings to non-dairy holdings, and within the non-dairy sector. Spread was least common from non-dairy holdings to dairy holdings. The number of secondary cases stratified by affiliation to the existing MAP voluntary surveillance programme is presented in Fig. 4 . Spread from programme member holdings to non-member holdings resulted in more secondary cases relative to all other scenarios. In contrast, no secondary cases were observed from non-member holdings to holdings affiliated to the programme.

Fig. 4

The median (solid black line), average (dashed black line) and the 95 % credibility interval (grey area) of the number of secondary infected holdings in the voluntary surveillance programme or non-programme member holdings that remain infected over time after an introduction of MAP infection into a single either member or non-member holding.

Surveillance performance

The estimated sensitivity of bulk milk ELISA testing that was used in the model is presented in Fig. 5 . The regression coefficient for -log2(within-holding prevalence) in this model was -2.03 (95 % confidence interval: -2.48, -1.59) with an intercept of 7.95 (95 % confidence interval: 6.12, 9.79) log odds of bulk milk test positivity.

Fig. 5

The solid black line represents the model predicted bulk-milk sensitivity from a logistic regression model of the proportion of animals that are positive in the holding. The grey area is the 95 % confidence interval around the model prediction. Data extracted from Nielsen and Toft (2014)Fig. 3.

When this bulk milk testing was applied to a scenario where MAP was introduced into a single holding, the performances of the various tested surveillance scenarios were evaluated and are presented in Fig. 6 . The outbreak died out within the 13-year model period in 60.1 % of scenarios and in those scenarios where the outbreak persisted, the proportion of detected outbreaks is summarised in Table 1 . The size of the outbreak at the time of detection and size of the missed outbreaks tended to decrease as the number of samples collected increased but did not differ between random holding-selection and targeted selection of holdings based on number of ingoing contacts due to movement of animals. These findings and the size of outbreaks missed by the surveillance are described in more detail in Appendix 1 - Supplementary material. The Cox regression analysis indicated that the time to detection was significantly different between surveillance methods. The hazard ratios (HR) for number of samples with 500 holdings as the reference group, were: 1.67 (95 % CI: 1.60–1.75) for 1000 holdings, 2.20 (95 % CI: 2.10–2.30) for 1500 holdings and 2.58 (95 % CI: 2.47–2.70) for 2000 holdings. The hazard ratio for holding selection by log2(indegree) was 0.94 (95 % CI: 0.91−0.97) and indegree targeting was 0.81 (95 % CI: 0.78−0.83) with random selection as the reference group. Results of this model indicate that the differences between holding-selection strategies was generally very small for efficacy of detecting point introductions over the first 13 years of an outbreak. In 95 % of outbreaks, ≤ 6 holdings were infected after 13 years and 50 % of outbreaks died out after 7.3 years.

Fig. 6

Simulated outbreak and surveillance progression given a seed of a single holding in the population infected at a within-holding prevalence of 2%, evenly distributed over calves, heifers and cows. The top graph indicates the median (solid black line), mean (dotted black line) low and high 95 % credibility interval of between-holding prevalence (grey area). The bottom graphs show the Cox proportional hazards model predicted cumulative hazard interpreted as the expected number of detections over time. The solid black line shows the mean number of detected holdings by random sampling, the solid red line shows the mean number of holdings detected when sampling was targeted by log2(indegree), the solid blue line shows the mean number of holdings detected when sampling was targeted by indegree. The corresponding coloured shaded areas are the upper and lower 95 % confidence intervals of the model estimates. Note that the median number of detected holdings over the 13-year period was 0 for all surveillance strategies, hence the cumulative hazard below 1.

Table 1

The percentage of persistent MAP outbreaks starting from a point introduction that were detected by surveillance systems with three alternative holding selection strategies over the study period. Simulated outbreaks and surveillance progression given a seed of a single holding in the population infected at a within-holding prevalence of 2%, evenly distributed over calves, heifers and cows. The outbreaks were considered to be persistent if holdings were infected on the last time step of the simulation. Outbreaks were persistent in 39.9 % of 342,024 model trajectories.

Number of samples collected per year	Selection of holdings targeted by indegree	Selection of holdings targeted by log2(indegree)	Random selection of holdings for surveillance
500	6.6 %	7.8 %	7.8 %
1000	10.3 %	12.4 %	12.7 %
1500	13.8 %	15.4 %	16.5 %
2000	15.7 %	17.9 %	18.7 %

When bulk milk testing was applied to a scenario where MAP was seeded to mimic a low prevalence endemic state of 0.2 % of holdings, the number of infected and detected holdings from bulk milk surveillance is presented in Fig. 7 . The Cox regression analysis indicated that the time to detection was significantly different between surveillance methods. The HRs for number of sampled holdings with 500 holdings as the reference group, were: 1.78 (95 % CI: 1.69–1.88) for 1000 holdings, 2.46 (95 % CI: 2.33–2.60) for 1500 holdings and 3.12 (95 % CI: 2.95–3.29) for 2000 holdings. The hazard ratio for holdings selection by log2(indegree) was 0.91 (95 % CI: 0.87−0.95) and 1.04 (95 % CI: 0.99–1.08; P = 0.114) for indegree targeting with random selection as the reference group. The proportion of detected trajectories with the endemic state of 0.2 % of holdings is summarised in Table 2 .

Fig. 7

Simulated outbreak and surveillance progression given a seed of 0.2 % of holdings in the population infected at a within-holding prevalence of 2%, evenly distributed over calves, heifers and cows. The top graph indicates the median (solid black line), low and high 95 % credibility interval of between-holding prevalence (grey area). The bottom graphs show the Cox proportional hazards model predicted cumulative hazard interpreted as the expected number of bulk-milk detections over time. The solid black line shows the mean number of detected holdings by random sampling, the solid red line shows the mean number of holdings detected when sampling was targeted by log2(indegree), the solid blue line shows the mean number of holdings detected when sampling was targeted by indegree. The corresponding coloured shaded areas are the upper and lower 95 % confidence intervals of the model estimates. The time where the cumulative hazard significantly exceeds 1 can be interpreted as the time to detection.

Table 2

The percentage of trajectories (12,000) detected by surveillance systems with three alternative holding selection strategies over the study period. Simulated disease spread and surveillance progression given a seed of 0.2 % of holdings in the population infected at a within-holding prevalence of 2%, evenly distributed over calves, heifers and cows. All trajectories had persistent MAP throughout the study period.

Number of samples collected per year	Selection of holdings targeted by indegree	Selection of holdings targeted by log2(indegree)	Random selection of holdings for surveillance
500	88.8 %	85.5 %	88.8 %
1000	97.7 %	97.0 %	97.1 %
1500	99.5 %	98.3 %	99.1 %
2000	99.7 %	98.7 %	99.3 %

For both Cox regression models, the residuals were flat over time indicating proportional hazards and a plot of the Cox-Snell residuals against the cumulative hazard followed a diagonal indicating a good model fit.

Discussion

This study illustrates how modelling can be used to plan and compare strategies for disease surveillance in populations where the disease is apparently absent, and the detection capacity of the diagnostic test varies by within-holding prevalence. The modelling framework included three components - the population demographics, the time-varying contacts between holdings, and disease testing. This allowed for estimation of disease spread and surveillance performance simultaneously and could be used to inform decisions about resource allocation and support the design of more cost-effective surveillance.

The model of surveillance in a low prevalence endemic state provides an assessment of how useful such a surveillance strategy would be if there existed a small fraction of holdings in the population that had not been detected by the previous systematic investigations or the ongoing clinical surveillance in Sweden. When holdings were seeded at the beginning of the study period, there was a short phase where the between-holding prevalence decreased. The infection died out in some holdings in this phase perhaps due to being isolated from others in the network as was illustrated in the within-holding model. It should also be noted that the number of holdings decreased considerably in Sweden during the study. The model settled on approximately 40 infected holdings which could then be detected by bulk milk surveillance during the 13-year study period. As more samples were collected, the surveillance strategies performed better by detecting infected holdings more frequently and faster. Even when the surveillance strategy only included 500 samples per year, the system most frequently detected at least one case within the 13-year period (Table 2). If 2000 samples were tested per year, then there was a 95 % probability of detecting at least 1 infected holding in the population within 3 years. This result illustrates the potential of using a very low sensitivity test on bulk milk to detect occurrence of disease in a large population that otherwise might not be useful to draw conclusions about a single holding.

In the scenario where MAP is absent from the population, a surveillance system should have the capacity to detect an introduction before the disease spreads to many holdings. The current modelling indicates that, for a slow spreading disease in the Swedish cattle population this capacity is difficult to assess over a 13-year period. Most frequently, the outbreaks died out or affected a very small number of holdings. The surveillance systems tested were also not sufficient to detect the remaining outbreaks that persisted. Even if 2000 samples were collected per year, the system only detected 0.2 positive holdings on average, and in most cases no holdings were detected. It is apparent that the prevalence of MAP must increase to levels that take longer than 13 years to develop after a single holding introduction before it is detectable using the current strategies. However, as the model of the endemic state indicated, MAP is readily detectable by bulk milk surveillance if the prevalence exceeds 0.2 % of the population.

When designing surveillance sampling, targeting higher risk units is almost universally accepted as a profitable strategy when the goal is to more quickly detect disease or demonstrate freedom in a population. In the current work, holdings were targeted randomly and by their history of introducing animals from other holdings. Perhaps counterintuitively, it was observed that the best strategy was to sample holdings at random or that risk-based sampling was equivalent to random targeting. Previously, targeting holdings with a high indegree has been shown to be an effective sampling strategy in the Swedish cattle population for the much faster spreading bovine coronavirus and bovine respiratory syncytial virus (Frössling et al., 2012). Under a targeted surveillance strategy, holdings with a high indegree were sampled as frequently as once per year. Whereas, in the random sampling scheme, these holdings were rarely sampled more than once over the 13-year study. In the case of a very slow spreading disease like MAP, the advantage to re-sampling high-risk holdings may be outweighed by the loss of coverage of the population. The log transformation of indegree reduced this resampling of holdings but still did not result in better surveillance performance. It is also plausible that because the spread of MAP in the current simulations was initialized randomly and the spread was very slow that this favoured random sampling. If spread was allowed to persist over longer than the 13-year study period the higher indegree holdings might be more likely to be infected and therefore more valuable to target in a surveillance.

The relationship between diagnostic sensitivity and the within-holding prevalence has been investigated for several diseases and specifically for MAP (Van Weering et al., 2007; Nielsen and Toft, 2014; Collins et al., 2017). From these studies, it is known that the ability of a bulk milk test for detection of antibodies will vary by the within-holding prevalence, and we therefore set up our simulations to account for this. Nielsen and Toft (2014) concluded that in practise, bulk milk testing for MAP has limited utility to make conclusions about the herd status of the disease. However, for testing a population of many holdings to demonstrate disease freedom and continuous detection of disease introductions, we show that although bulk milk testing often has a low holding sensitivity, it could still be a useful tool for surveillance of MAP. In the current work, false positive test results are not considered but would need to be managed if a testing scheme were implemented. Bulk milk testing was investigated in the current study for the practical reason of access to these samples as part of various other surveillance programmes currently applied in Sweden. The use of bulk milk exclusively for surveillance of MAP would result in coverage that would not directly include non-lactating cattle (calves and dry cows) or non-dairy herds but rely on the potential spread to those animals and herds that are under surveillance before detection is possible. However, if these samples could be used to also monitor MAP in Sweden, the additional cost of sample collection could also be avoided. In the current study, the sensitivity of bulk milk testing (SeBM) was derived by a two-step approach using estimates of sensitivity of the individual (SeI) from Nielsen et al. (2013) and the age distribution of Swedish cattle to determine apparent within herd prevalence and subsequently evidence from Nielsen and Toft (2014) to estimate SeBM for each holding at the time of testing. The case definition in Nielsen et al. (2013) was based on testing and not the absolute state of infection as is available in the current model and this could have resulted in inflated estimates of SeI and subsequently SeBM due to misclassification bias.

Given the population structure and contact network of the Swedish cattle population, our model results suggest that potential spread of MAP would be slow and partly limited by barriers in the contacts between sectors in the cattle industry between non-dairy, dairy and holdings that voluntarily restrict trade with holdings outside an existing voluntary surveillance programme. MAP has never been detected in Swedish dairy cattle and no dairy cattle have been imported since Sweden joined the EU in 1995. To ensure that introductions of MAP remains detectable in the different subpopulations, more active surveillance activities targeting both non-dairy and dairy cattle may be preferable. However, the lack of overlap and slow between-holding spread is, of course, also a beneficial protective factor that decreases the probability of any future introductions of MAP spreading to the dairy population.

The initialisation of the model by seeding a holding or multiple holdings with the rounded up 2% within-holding prevalence across all age categories was chosen to try to start an outbreak without it immediately dying out. This was not meant to mimic a plausible introduction or an endemic state only a means to seed the model with infected animals. One might argue that seeding with fewer animals or only in certain age categories to simulate importing a smaller number of infected animals to the population would be an interesting scenario to model and test surveillance against. However, this would result in many more trajectories of the model where the disease simply dies out with no chance to determine the performance of the subsequent surveillance. The seeding of disease in the model was also done at the beginning of the study period in order to maximize the available study time due to the limited access to 13 years of animal movements. Seasonality of animal movements could therefore have introduced a bias in the results because the seeding always occurred at the same time of year.

The within-holding model for spread of MAP used in the current work was adopted from a previous study that was parameterised to fit the study of the spread of MAP over 10 years in 3 US dairy herds. Since spread of any disease is affected by herd management, animal flow, herd size, etc., the parameters in the model are both affected by the biology of MAP and management factors in these US dairy herds. Due to these potential differences in production systems, we must consider that the spread of MAP in a Swedish herd could be different from what this model predicts and therefore augment the uncertainty of the final model estimates. Others have also modelled the spread of MAP including a highly detailed mechanistic within-herd model (Kirkeby et al., 2016) and other compartments models that include a between herd spread (Beaunée et al., 2015; Magombedze et al., 2013) as well as scenario tree models of surveillance efficacy (Meyer et al., 2018; Frössling et al., 2013). These models and approaches were considered for the current work but the structures could not readily be used in the SimInf model framework in order to the include within- and between-herd animal movement events as was done in the current work to precisely reflect the Swedish cattle industry. The inclusion of these movement events allowed the model to reflect the changing population structure over time such as reduction in the number of herds and increasing herd size over time as well as seasonal changes in animal distributions and the network connections in the industry. These features of population dynamics and network connections are understood to play an important role in disease spread and could in this way be included to mimic the Swedish cattle industry without the need to estimate a synthetic model of this complex network.

The use of a compartment model as a tool, not only to model the change in states of animals in the holdings, but also to keep track of the states of testing holdings in the model proved to be very useful. It allowed the surveillance to be modelled in real-time alongside the disease spread and to be pre-scheduled, similar to the way surveillance is designed in practise. Summarising the data from simulations was challenging since the surveillance system sensitivity and timeliness were coupled. To simply summarize the time to detection in those outbreaks that were detected would result in bias due to differences in sensitivity, and survival analysis which could account for censoring was thus required. This type of adjustment was not done for the size of outbreak at the time of detection which may have been the reason that differences were not detected. However, the probability of detection of an outbreak is also associated with outbreak size, which also makes this parameter inherently difficult to interpret. Currently, single samples were tested in each holding, but the approach could also be generalised to multiple samples per holding in the case of modelling individual-animal tests. Also, the findings from the surveillance were only recorded as observations from the model. In future implementations we see the potential to model complex disease spread and surveillance systems incorporating surveillance results as feedback to the disease spread model. This could allow for testing intervention strategies in the model to reduce prevalence of disease in the population over time, not just detect it.

Conclusion

This study describes how simulation of surveillance in a disease spread model can be used to assess several surveillance-performance parameters, and compare different surveillance strategies. The results indicate that, although the diagnostic sensitivity is varying, testing for presence of MAP antibodies in bulk milk can be a useful in surveillance efforts performed to demonstrate disease freedom. The results further suggest that the spread of MAP in the Swedish network of cattle holdings can be expected to be slow and partly limited to certain subpopulations based on e.g. industry sector. The use of a disease spread modelling approach to simulate surveillance allowed for the evaluation of surveillance after adjusting for complex changing population structure and concurrent disease spread.