Infection seeding

We seeded the epidemic with five infected cows in a mid-size herd in Texas, on the week beginning December 18th 2023, based on phylogenetic analyses⁶. For the stochastic realizations, we also seeded 9 additional herds in accordance with the nine early outbreaks detailed in Caserta et al.³. The herd size, number of infected cattle, and date of seeding is consistent with the data presented in that manuscript.

Epidemiological dynamics

We construct a stochastic metapopulation SEIR model³⁶ with 35,974 individual herds of varying population size, informed by the 2022 US Agricultural Census¹¹. Each herd’s infection dynamics are the stochastic equivalent of the following set of ordinary differential equations (ODEs):

Here, S_i, E_i, I_i, and R_i are the number of susceptible, exposed, infected and recovered cows in herd i. N_i is the total population of herd i. β, σ, and γ are the transmission, incubation, and recovery rates respectively. α is a model parameter between 0 and 1 controlling the rate of transmission between herds in the same state. I_−i and N_−i are the total number of infected cattle, and the total number of all cattle, in the US state herd i resides in, not including the cattle in herd i itself. Early epidemiological surveys of farms reporting outbreaks found that transmission routes existed between herds in the same state through the shared use of equipment, staff, or the movements of wild birds³⁷, which we capture here in the model. We assume no such forms of transmission can occur between herds in different US states.

The stochastic analogue of the above ODEs, is that we calculate the number of cattle progressing between epidemiological compartments via binomial distributions, for each time step dt as:

Here \({n}_{XY}^{i}\) is the number of cattle moved from compartment X to Y (for general X and Y), in herd i, in a time step of size dt.

After all cattle movements between epidemiological compartments is concluded, we calculate for each herd that has yet to report an outbreak, whether or not it will report an outbreak in that time step. It reports an outbreak with probability \({P}_{i}^{\,{\mbox{outbreak}}\,}=1-{e}^{-{\phi }_{i}}\), where ϕ_i is

and A^asc is a model parameter that we fit. The bracketed term to the left of A^asc in Eq. (3) is shown in the heatmap of Fig. 3A. This functional form was developed in consultation with veterinarians based on their experiences of at what stage of pathogen spread they are typically consulted. While US states undoubtedly vary in their detection capabilities, there is insufficient outbreak data to fit unique A^asc values for each state. Assuming one national A^asc parameter allows us to identify which states that have reported 0 outbreaks to date are driven mostly by under-reporting (Fig. 2B).

Movement of cattle between herds

After calculating the movement between epidemiological compartments and any reporting of outbreaks, we then calculate the movement of cattle between herds. As detailed in Supplementary Material Section 2.4, we infer from the USAMM the probability, \({P}_{k}^{\,{\mbox{export}}\,}\), for each US state, k, that a herd within that state will export cattle each week. We assume the same probability for every herd in the state. We also calculate the proportion of cows in the origin herd that will be exported—\({P}_{k}^{\,{\mbox{export size}}\,}\) from the USAMM export simulations, which include cohort size and size of origin herd. We also calculate the probabilities of, should an export of cattle occur, which US state they will be exported to. This is parameterized by a movement matrix M, where element M_k,l denotes the probability that an export from state k will go to state l. This matrix describes the patterns of interstate movement, and the diagonal represents the probability of an export remaining within the same state. The exact matrix is provided as Supplementary Data. Once the destination state is determined, we randomly allocate which herd in the destination state the cattle will be exported to, scaled by the population size of the respective herds, to preserve herd sizes. Once an origin herd, i, and destination herd, j, are assigned, we draw the number of cattle to be exported as

where k is the US state that origin herd i resides in. Lastly, before moving cattle between the respective compartments of herds i and j, we simulate the border testing mandate. If the model date is after April 29th 2024, we draw a random variable, X from a hypergeometric distribution:

Here the three parameters of the above hypergeometric are, the number of success items in the population, the number of failure items in the population, and the number of samples taken without replacement from the population. X is the number of infected cattle drawn. If X = 0, then no infected cattle are detected, and the export takes place. Note, a positive test prevents the export, but does not immediately register as a reported outbreak. All probabilities and a full logic flow diagram are presented in Supplementary Material Section 2. U.S. state boundaries were obtained using the maps package in R (via map_data(“state”)) and visualized with ggplot2.

cowflu package

To efficiently simulate the above probabilistic model, we produced a custom R package, cowflu³⁸, which allows simulating and fitting the model via the dust2 package²² in R, while the model itself is written in C++. Documentation on the use of the package and worked vignettes can be found on our github repo: https://github.com/mrc-ide/cowflu. The package is flexible to being applied to any SEIR metapopulation model with custom probabilities of movement between sub-populations, subject to user-defined movement matrices.

Model fitting

Five of the above model parameters—β, α, σ, γ, and A^asc, are fit via particle Markov Chain Monte Carlo²⁴ methods. We assign weakly-informative prior distributions, informed by early studies associated with the current outbreak³⁹. We fit the model simulated values of date of first outbreak detection (as seen in Fig. 2A) to the real world data equivalent, via a likelihood function detailed in Supplementary Material section 2.5. We ran the pMCMC simulations across 16 chains of 40,000 iterations each. Model convergence statistics are presented in Supplementary Material section 2.5.

Table 2 shows the priors and posteriors for all model parameters. Note that we fit \(\frac{\beta }{\gamma }\) instead of β due to observed correlation between β and γ, so as to improve chain mixing.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.