PRRS Outbreak Prediction via Deep Switching Auto-Regressive Factorization Modeling

The Porcine Reproductive and Respiratory Syndrome (PRRS) is arguably the most challenging and costly viral infectious disease in the US swine industry [1]. It is the high cost of disease control practices such as testing and vaccination as well as the financial damage caused by an outbreak that highlights the need to develop predictive models that can help to identify farms at high risk of infection to support risk-based, more cost-effective, target interventions. Such a framework will allow for more efficient testing, vaccination and outbreak prevention. Due to the high level of specialization in the swine industry, vast amount of data has been collected by the swine production systems. However, they have not yet been exploited enough due to difficulty of data access, integration and analyses (i.e., data are not consistently gathered, are non-standardized which makes their integration difficult, and are usually scattered across stakeholders). Examples of these data include diagnostic information, including the number of infected or dead animals, animal movements between farms or production data, which can give insights regarding farm health status and its contact network. Usually, these data do not satisfy the granularity required for learning an advanced predictive model (e.g. diagnostic samples are only taken once or twice per month per farm). However, using the real-world data, we can simulate epidemics to produce fine-grained time series data to analyze it further with an advanced novel prediction method based on a generative and variational inference model.

The direct contact with an infected pig is the main pathway for PRRS virus transmission. The network of between-farm movements, shown in Fig. 1a together with the intra-farm (local) pig-level contact pattern based on Susceptible-Exposed-Infectious-Removed (SEIR) epidemic model (see Fig. 1b) allows for constructing a system-level (global) pig-level disease transmission contact network, [2]. Network-based SIR or SIR-extended epidemic models have been extensively studied in the literature [3]. In [4], Newman studied a network-based SIR epidemic model where infection is transmitted through a random network of contacts between individuals. The disease transmission contact network is a probabilistic graph that can be sampled to generate virtual contact. Using this, we generate fine-grained spatio-temporal time-series data based on statistics of real-world data.

The spatio-temporal data is often considered to have a high level of correlation between spatial dimensions, and, therefore, they can be assumed to be governed by a smaller number of underlying components. For modeling the temporal dynamic of the time series including the number of infected, dead, or recovered pigs, we employ a non-linear vector auto-regressive latent model inspired by the work of [5]. Our spatio-temporal time-series data is first factorized into temporal weights and spatial factors. The temporal weights are modeled using a non-linear auto-regressive model parameterized by neural networks governed with a Markovian chain of discrete switches to capture higher-order multimodal latent dependencies, [6, 7, 8, 9, 10].

Based on the rich database of an extensive anonymous swine production system located in the Midwest of the United States, we have access to farm-level pig shipment data, and PRRSV testing results [11]. From 2006 to 2021, there have been over 260,000 movement records to or from farms within farm farm this production system. For each movement entry, the data includes the source and destination information, the number of transported pigs, and the date of the movement. Based on the farm-level shipment data we generate a farm-level movement network for the entire production system. Furthermore, the frequent PRRSV testing in each farm gives insight into how the virus is transmitted, e.g., what is the virus’s transmission rate, incubation time, etc. Using the SEIR model we can produce a pig-level contact network, [2]. The combination of this two network built on the statistics extracted from our real data results in an intricate contact network by which we can simulate complex time series data showing the number of infected, dead, or recovered pigs in each farm in the production system.

To create the farm-level movement network we build a probabilistic graph, i.e., a graph in which the existence of edges is uncertain with some probability. A node in the graph represents a farm, while the weight of an edge is proportional to the shipment rate ( probability of shipment) between farms. Fig. 1a shows a part of the shipment network of over 300 farms for our real data. The edge thickness is representative of the shipment rate. For local (intra-farm) pig-level contact network in each farm, we consider a basic random graph model based on Erdös—Rényi model [12, 2], that produces pig contact graphs with an edge probability of $0.5$ between any pair of pigs. The global (inter-farm) pig-level contact network is constructed when we sample a random generalization of a between-farm shipment over time (each day) and as a result, we create pig contact between farms (see Fig. 1b, top).

The simulation is formed on the network-based SEIR epidemic model for PRRS. In a network-based model, we consider a graph in which nodes represent individual pigs, and edges indicate direct or indirect contacts between pigs, which are considered infection pathways of PRRS. Each animal can be in one of the four states, Susceptible (S), Exposed (E), Infected (I), or Recovered (R) as the result of the epidemic progression. The state-transition diagram between these states is shown in figure. 1b, Bottom. In the generated swine pig-level network, a PRRS outbreak is introduced by randomly selecting a pig farm, and infecting an arbitrary random number of pigs. We collect several time snapshots representing the progression dynamics of the disease spread, such as the number of infected pigs in each farm over time. The healthy pigs which are free from PRRS virus infection are classified as Susceptibles. If such a healthy pig comes into contact with infected pigs containing the virus, it may get infected at the rate $\beta Y_{i}(t)$, where $Y_{i}(t)$ is the number of infected neighbors of node $i$ at time $t$. If the transmission of pathogen occurs, a healthy pig enters into the Exposed group where it stays for the duration of the incubation period. On average, this period is denoted by $\nicefrac{{1}}{{\sigma}}$. Once it shows symptoms, it moves into the Infected group. It stays there for an average time of $\nicefrac{{1}}{{\gamma}}$ before it is recovered. We choose the parameters values $\beta=0.087,\sigma=7,\gamma=6.5$ based on [13, 14]. This dynamic produces a spatio-temporal time series from all farms over time that can be used for predication modeling in the next section (Sec. 3).

Our spatio-temporal data indicates the number of pigs categorized within a particular stage, e.g., infected, recovered, etc., in every time instance in each farm. We denote this data as the matrix $X\in\mathbb{R}^{T\times D}$, where $T$ is the number of time points and $D$ the number of spatial locations, e.g., the number of farms. Building on previous work by [5], our assumption is that $X$ can be decomposed into a weighted summation of $K\ll D$ factors over time as:

$$X\approx[w_{1},\cdots,w_{T}]^{\top}[f_{1};\cdots;f_{K}]=W^{\top}F,$$

(1)

where $f_{k}\in\mathbb{R}^{D}$ is the $k^{\text{th}}$ spatial factor and $w_{t}\in\mathbb{R}^{K}$ is the weight vector at time $t$. Our intuition for adopting this model for some pig-specific collected measurements in $D$ farm over $T$ time points is that there are $K\ll D$ underlying factors using which we can approximate the overall dynamics of the disease propagation in the data. Our main goal is to predict the future outbreak behavior given the past time samples in each farm.

We assume that the weights, $W=\{w_{t}\}_{t=1}^{T}$, are generated according to a set of temporal lags, $\ell$, through a deep probabilistic switching auto-regressive model. These weights are furthermore governed by a Markovian chain of discrete latent states, $\mathcal{S}=\{s_{t}\}_{t=1}^{T}$ as follows: $w_{t}\sim p(w_{t}|w_{t-\ell},s_{t}),\;s_{t}\sim p(s_{t}|s_{t-1})$. In addition, we assume that spatial factors, $F=\{f_{k}\}_{k=1}^{K}$, are controlled by a shared low dimensional latent variable, $z$, as follows: $f_{1:K}\sim p(F|z),\;z\sim p(z)$. Fig. 2 shows the relation among above-mentioned random variables in a probabilistic graphical model diagram form.

We train the model using stochastic variational methods [15, 16, 17, 18] by approximating the posterior $p_{\theta}(\mathcal{S},W,z,F|X)$ using a variational distribution $q_{\phi}(\mathcal{S},W,z,F)$, and by maximizing a lower bound (known as ELBO) $\mathcal{L}(\theta,\phi)\leq\log p_{\theta}(X)$:

	$\displaystyle\mathcal{L}(\theta,\phi)$	$\displaystyle=\mathbb{E}_{q_{\phi}(\mathcal{S},W,z,F)}\left[\log\frac{p_{\theta}(X,\mathcal{S},W,z,F)}{q_{\phi}(\mathcal{S},W,z,F)}\right]$		(2)
		$\displaystyle=\log p_{\theta}(X)-\text{KL}(q_{\phi}(\mathcal{S},W,z,F)\,||\,p_{\theta}(\mathcal{S},W,z,F|X)).$

By maximizing the bound with respect to the parameters $\theta$, we learn the generative distribution over datasets $p_{\theta}(X)$, and by maximizing the bound over the parameters $\phi$, we do Bayesian inference by approximating the distribution $q_{\phi}(\mathcal{S},W,z,F)\simeq p_{\theta}(\mathcal{S},W,z,F|X)$ over latent variables for each data point. According to the graphical model in Fig. 2, the joint distribution of observations and latents will be:

	$\displaystyle p_{\theta}(X,\mathcal{S},\mathcal{Z})=p(F|$	$\displaystyle z)p(z)\prod_{n=1}^{N}p(X_{n}|W_{n},F)p(w_{n,-\ell})p(s_{n,0})$
		$\displaystyle\prod_{t=1}^{T}p(s_{n,t}|s_{n,t-1})p(w_{n,t}|w_{n,t-\ell},s_{n,t}),$		(3)

where $\mathcal{Z}=\{W,z,F\}$). Furthermore, we assume a fully factorized variational distribution for the latent variables posterior as:

$\displaystyle q(\mathcal{S},\mathcal{Z})=q(F)q(z)\prod_{n=1}^{N}q(w_{n,-\ell})q(s_{n,0})\prod_{t=1}^{T}q(s_{n,t})q(w_{n,t}).$

(4)

We used time series of epidemic progression from over 300 farms simulated for 700 time points. We kept last $20\%$ of the time series as the test set. We then performed a short-term prediction tasks by adopting a rolling prediction scheme reported in [21]. For short-term prediction, the next time point is predicted on the test set using the generative model and spatial factors learned on the train set. We reported the test set normalized root-mean-square error (NRMSE%), which is related to the expected negative test-set log-likelihood for the case of Gaussian distributions, and it is used for evaluating the predictive generative models. We obtained NRMSE of $2.5\%$ averaged over all the farms. Fig. 3 shows the number of infected pigs over time for nine selected farms. In this figure, we illustrate the actual number of infected pigs in the simulated data with a solid green curve, the mean estimate of the predictive model with a dashed purple curve, and the standard deviation of the data with a shaded red error bar. Each row in the figure represents a group of relatively highly connected farms in terms of the frequency of the pig shipments. Note that each group show a relatively strong correlation regarding the outbreak progression and the predictive model was able to capture this.

One observation regarding the performance of the model in situations when we do not have a curve with smooth behavior is the fact that when the increase and decrease of the number of infected pigs is abrupt, the model tends to converge to a middle point between the current value and the future one. This issue can be addressed using a model selection choice. Specifically, we can control the number of factors $K$ in order to increase the flexibility of the model in capturing different levels of curve smoothness.

The PRRS outbreak cause an economic loss of over $664 million annually [1], which can be significantly mitigated by early detection and risk-based intervention practices. Direct contact is the main disease transmission pathway. Therefore, the pig contact network provides a substantial basis to develop an outbreak prediction framework. We create a system-wide pig contact network by combining the SEIR epidemic model based on intra-farm infection transmission parameters, and inter-farm pig shipment network. We presented a hierarchical factorized deep generative model of our spatio-temporal data that can capture the underlying dynamics of the disease spread with the aim to predict the number of infected pigs in all farms. Our result demonstrates the ability of the model in forecasting the virus spread progression with an average one-day prediction error of NRMSE = $2.5\%$. A potential future direction is to incorporate the per farm disease transmission parameters for the SEIR model to represent variations in the disease control implementation. Additionally, indirect disease transmission pathways, such as airborne, can be included in the framework for a more comprehensive analysis.