Complex responses to movement-based disease control: when livestock trading helps

Livestock disease controls are often linked to movements between farms, for example, via quarantine and pre- or post-movement testing. Designing effective controls, therefore, benefits from accurate assessment of herd-to-herd transmission. Household models of human infections make use of R*, the number of groups infected by an initial infected group, which is a metapopulation level analogue of the basic reproduction number R0 that provides a better characterization of disease spread in a metapopulation. However, existing approaches to calculate R* do not account for individual movements between locations which means we lack suitable tools for livestock systems. We address this gap using next-generation matrix approaches to capture movements explicitly and introduce novel tools to calculate R* in any populations coupled by individual movements. We show that depletion of infectives in the source group, which hastens its recovery, is a phenomenon with important implications for design and efficacy of movement-based controls. Underpinning our results is the observation that R* peaks at intermediate livestock movement rates. Consequently, under movement-based controls, infection could be controlled at high movement rates but persist at intermediate rates. Thus, once control schemes are present in a livestock system, a reduction in movements can counterintuitively lead to increased disease prevalence. We illustrate our results using four important livestock diseases (bovine viral diarrhoea, bovine herpes virus, Johne's disease and Escherichia coli O157) that each persist across different movement rate ranges with the consequence that a change in livestock movements could help control one disease, but exacerbate another.


General SIS model
Note that we exclude events which have no effect on the state, e.g. all deaths are coupled to replacement by new susceptibles, so the state change for mortality of S is S → S. In later models, we do not specify every single movement combination, and instead just note the general case X → Y at rate κP Y X.
In all cases where the SIS model was used, mortality µ = 1/3 [2], while values used for β and γ are specified in the relevant figures.

Escherichia coli O157
We used an SLHS model, where susceptibles S become either super shedders H (high), or regular shedders L (low), with probabilities p and 1 − p respectively. Super shedders are η times more infectious than regular infectives.
and where the parameters are as specified in Table 1. We took R 0 = 1.5 [3], and assumed a mean infectious period of 4 weeks. The stochastic event table is shown in Table 2.
We calculate R * for this model numerically in Section 3.1 of the main manuscript using the simulation approach discussed in Section 2, but we also calculate R 0 using the NGM method, to show that R 0 does not depend on the probability of super shedding p.
Using the NGM method, the Jacobian matrix of the linearised subsystem is and the transition matrix and so the NGM is Therefore we choose η = (1/p − 1) 2 to ensure the correct transmissibility, and c = 1/(p − 1) to normalise transmission, and then R 0 = β/(µ + γ) is independent of p.

Bovine herpes virus (BHV)
We used an SILI model, where susceptibles S become infectives I, then recover to a latent state L, from which they may later relapse to I again.
and where the parameters are as specified in Table 1. The stochastic event table is shown in Table 2.

Bovine viral diarrhoea virus (BVDV)
We used an SIPR model, where individuals start as susceptibles S, then become infectives I, and recover to resistant R. Newly born individuals can become persistent shedders P with vertical transmission probability p v .
and where the parameters are as specified in Table 1. The stochastic event table is shown in Table 2.

Mycobacterium avium ssp paratuberculosis
This is the pathogen responsible for Johne's disease, and is also known as Map or ParaTB. We used an SLHC model, where susceptibles S become low shedders L, then progress to high shedders H, and finally progress to clinically infected C. Newly born individuals become low shedders L with vertical transmission probability p v .

Simulation methods
From an ease of computation perspective, there are several approaches to calculating R * , depending on model complexity and disease dynamics. These methods may involve calculation of T inf and P pos , or by counting the number of infectives leaving the primary herd. The Master Equation [4] is a set of differential equations that describe the time evolution of the probability distribution of all possible states in a stochastic model. The stochastic simulation algorithm [5] (also known as the Gillespie algorithm) provides samples which are identically distributed to the solution of the Master Equation.
The Master Equation can be used to calculate R * exactly, by determining the elements T j inf and P j pos (i), the time infection persists following introduction of an individual of infectious class j, and the average prevalence of infectious class i during that period, then using these to populate the NGM K as described in Section 2 of the main manuscript. All figures for the SIS model were generated using the exact solution from the Master Equation.
However, the Master Equation rises in complexity with the number of possible states, and is not appropriate for models with multiple disease states and large herd size N . For example, with N = 50, the SIS model has 51 possible states, while the SLHC model for ParaTB has 22,100 possible states, with a correspondingly large Jacobian matrix, which can only be stored on most desktop computers as a sparse matrix.
For complex models it is more straightforward and biologically intuitive to use the stochastic simulation algorithm, proceeding straight to the NGM formulation and populating K directly. Each entry K ij is obtained by seeding the herd with one individual corresponding to type j, and observing the number of individuals of type i leaving the herd through movement before the disease becomes locally extinct. However, a large number of realisations of the stochastic model may be required to obtain a good approximation. All Figures in Section 4 of the main manuscript and Fig. 4a were calculated using this method.

Calculating the equilibrium proportion of infected herds
For the SIS model, iterative methods were used to determine the equilibrium proportion of infected herds, using the property that at endemic equilibrium, effective R * = 1 (i.e. at equilibrium, each infected herd infects on average exactly one other infected herd during its infectious lifetime). The equilibrium proportion may therefore be calculated by proposing a value x, which is the mean prevalence of infectives within the metapopulation P I , then increasing x if R * > 1, or decreasing x if R * < 1. Since 0 ≤ x ≤ 1, this may be effectively done using a binary search, until effective R * is sufficiently close to 1, and then taking the equilibrium proportion to be x/P pos .
For the livestock disease models, R * and equilibrium proportion were calculated via stochastic simu-lation, using a metapopulation model with n = 100 herds for the latter.