Bovine viral diarrhea (BVD) is a disease in cattle with complex transmission dynamics that causes substantial economic losses and affects animal welfare. The infection can be transient or persistent. The mostly asymptomatic persistently infected hosts are the main source for transmission of the virus. This characteristic makes it difficult to control the spreading of BVD. We develop a deterministic compartmental model for the spreading dynamics of BVD within a herd and derive the *basic reproduction number*. This epidemiological quantity indicates that identification and removal of persistently infected animals is a successful control strategy if the transmission rate of transiently infected animals is small. Removing persistently infected animals from the herd at birth results in recurrent outbreaks with decreasing peak prevalence. We propose a stochastic version of the compartmental model that includes stochasticity in the transmission parameters. This stochasticity leads to sustained oscillations in cases where the deterministic model predicts oscillations with decreasing amplitude. The results provide useful information for the design of control strategies.

Dynamical systems are often described by deterministic mathematical models, where the state of the system is determined by the initial conditions. Many real-world processes, however, include an element of probability or randomness. Even so, deterministic models might still be able to reproduce the main trend as a mean-field approximation but fail to capture the spectrum of possible dynamical scenarios of individual realizations. In addition, stochastic input such as noise can trigger the emergence of hidden dynamical features with surprising effects such as stochastic resonance, coherence resonance, or other noise-induced changes of dynamical behavior. Here, we present the example of a cattle disease that is realized as an extended susceptible-infected-recovered model. To explore the impact of stochasticity on the temporal behavior of the dynamics, we consider a stochastic transmission coefficient and systematically investigate the interplay between parameter noise and the intrinsic time scales of the underlying deterministic system.

## I. INTRODUCTION

Bovine viral diarrhea (BVD) is a viral disease that affects cattle and has a significant negative economic impact on the global livestock industry.^{1} BVD has a complicated pathogenesis that includes both transient (temporary) and persistent (life-long) infections. The spread of the bovine viral diarrhea virus (BVDV) occurs via both horizontal (contact between animals) and vertical (during certain stages of gestation) transmission.^{2} Acute infection in non-pregnant and non-immune cattle leads to a transient disease with complete recovery within 3 weeks.^{3} Clinical signs include fever, loss of appetite, mucosal lesions, and diarrhea^{4} with a very low associated mortality rate. The acute infection with BVDV induces a life-long protective immunity.^{5,6}

Vertical transmission, i.e., transmission from the mother to the fetus during pregnancy, is a complex process, which depends on the age of the fetus. Fetal infection in the period between around day 30 and day 120 can produce calves that remain persistently viremic for life.^{7} Abortions and teratogenic effects can result from infection during approximately the first 150 days.^{8–10} Calves infected during the last trimester have an active immune response.^{11} The full duration of bovine pregnancy is roughly 280 days.

Persistently infected (PI) animals lack an active immune response to the pathogen and excrete BVDV throughout their lives, and they are the most important sources of infection for BVDV. Possible symptoms include recurrent intestinal and pulmonary symptoms, neurological disorders, and growth retardation.^{12,13} Transiently infected (TI) animals are considerably less important as the source of infection,^{14} as they shed the virus in smaller quantities. PI cows may also give birth to persistently infected calves; however, fertility is reduced.^{15,16} The susceptibility of PI animals to mucosal disease results in high mortality,^{13} where mucosal disease, which is characterized by a mortality of almost 100%, develops only in PI animals. Clinical signs include anorexia and erosion of the intestinal tract and death follows approximately one week after the onset of symptoms.^{17} The pathogenesis of mucosal disease is not yet fully understood and is the subject of ongoing research.^{18,19}

There are various strategies that exist for controlling BVD. For a long time, control methods were limited to vaccination practices. This is a relatively inexpensive method but a successful strategy based solely on vaccination has never been reported.^{20} With growing knowledge about the pathogenesis of BVD and the development of diagnostic tests, PI animals have become the main target to control the spread of BVD and thus limit the associated economic impact.

Several models have been developed in order to study BVD within cattle populations. Some of them focus primarily on the estimation of economic loss or on various control measures,^{21–25} whereas others focus on spreading dynamics.^{26–29} Most models for BVD are discrete time stochastic models. One such model employs an agent-based approach,^{30} which allows the introduction of individual heterogeneities and complex network interactions. It is, however, difficult to derive analytical results from these agent-based models. Few authors have developed compartmental models with continuous time for the spreading dynamics of BVD.^{31–33} Innocent *et al.* describe a compartmental model and find broad agreement with a stochastic discrete model for large herd sizes.^{31} Basset developed a compartmental model formulated as integrodifferential equations.^{32}

In this paper, we investigate the stochastic effects of our BVD model and the impact they have on the spreading behavior of BVD. First, we present a deterministic compartmental model with continuous time that is based on a model suggested by Cherry *et al.*^{33} We identify steady-state solutions (equilibria) and analyze their stability in the context of a next-generation matrix. This enables the derivation of an insightful epidemiological quantity that characterizes the behavior of the spreading dynamics: the *basic reproduction number*, which quantifies the impact of an infected individual in terms of the number of expected secondary infections. Subsequently, we introduce a stochastic transmission coefficient and study its effect on the spreading dynamics.

## II. A DETERMINISTIC COMPARTMENTAL MODEL FOR THE SPREADING DYNAMICS OF BVD

### A. Model development

The model, schematically shown in Fig. 1, is now described step-by-step. The unit of the compartment variables is $hosts/km2$. There are six compartments considered in this model as follows:

$S$ is the fraction of susceptible animals. This compartment comprises three constant subgroups: (1) non-pregnant animals with fraction $p1$, (2) animals pregnant 1–150 days with fraction $p2$, and (3) animals pregnant 151–280 days with fraction $p3$ and $p1+p2+p3=1$.

$I$ denotes the fraction of transiently infected (TI) animals.

$P$ represents the fraction of persistently infected (PI) animals.

$R1$ describes the fraction of recovered, non-pregnant animals.

$R2$ denotes the fraction of recovered animals that were pregnant 1–150 days at the time of infection. After birth to a calf, they return to the $R1$ compartment.

$R3$ describes the fraction of recovered animals that were pregnant 151–280 days at the time of infection. After birth to a calf, they return to the $R1$ compartment.

Including the above-mentioned fractions $p1$, $p2$, $p3$, there are a total of nine parameters used in the model as summarized in Table I and discussed next. This discussion will finally lead to Eq. (3) below.

Parameter . | Value . | Definition . |
---|---|---|

β_{p} | 2.01/day | Transmission coefficient pertaining to |

PI animals | ||

β_{I} | 0.134/day | Transmission coefficient pertaining to |

TI animals | ||

μ | 0.001/day | Host birth/removal rate |

γ | 0.057/day | Recovery rate |

p_{1} | 0.45 | Probability that host is not pregnant |

when first exposed | ||

p_{2} | 0.3 | Probability that host is pregnant |

1–150 days when first exposed | ||

p_{3} | 0.25 | Probability that host is pregnant |

151–280 days when first exposed | ||

ϕ_{2} | 0.005/day | Reciprocal of the average time spent |

carrying an infected fetus | ||

ϕ_{3} | 0.021/day | Reciprocal of the average time spent |

carrying an actively immune fetus | ||

θ | 0.2 | Proportion of infected fetuses which |

survive to enter the herd | ||

a | 0.0008/day | Reduction in birth rate of PI animals |

b | 0.002/day | Additional mortality of PI animals |

Parameter . | Value . | Definition . |
---|---|---|

β_{p} | 2.01/day | Transmission coefficient pertaining to |

PI animals | ||

β_{I} | 0.134/day | Transmission coefficient pertaining to |

TI animals | ||

μ | 0.001/day | Host birth/removal rate |

γ | 0.057/day | Recovery rate |

p_{1} | 0.45 | Probability that host is not pregnant |

when first exposed | ||

p_{2} | 0.3 | Probability that host is pregnant |

1–150 days when first exposed | ||

p_{3} | 0.25 | Probability that host is pregnant |

151–280 days when first exposed | ||

ϕ_{2} | 0.005/day | Reciprocal of the average time spent |

carrying an infected fetus | ||

ϕ_{3} | 0.021/day | Reciprocal of the average time spent |

carrying an actively immune fetus | ||

θ | 0.2 | Proportion of infected fetuses which |

survive to enter the herd | ||

a | 0.0008/day | Reduction in birth rate of PI animals |

b | 0.002/day | Additional mortality of PI animals |

Animals in compartment $S$ move to the $I$ compartment due to interactions with both TI and PI animals; however, we assume this is independent of their pregnancy status. The transmission rate consists of two bi-linear terms corresponding to the infection rates from transiently and persistently infected animals giving the term $(\beta II+\beta PP)S$.

TI animals move to one of $R1$, $R2$, or $R3$ depending on their pregnancy status at the time of infection, i.e., whether they were part of the $p1S$, $p2S$, or $p3S$ susceptible subgroups, respectively, giving rise to the corresponding $\gamma p1I$, $\gamma p2I$, and $\gamma p3I$ terms.

Recovered animals in $R2$ give birth with rate $\varphi 2$ and, therefore, move into the non-pregnant recovered compartment $R1$. Likewise, recovered animals in $R3$ give birth with rate $\varphi 3$ and therefore move into the non-pregnant recovered compartment $R1$. Additionally, births from the $R2$ compartment become persistently infected and move into the $P$ compartment; however, due to infection, not all calves survive and the number of births $\varphi 2R2$ is reduced by the factor $\theta $, giving the term $\theta \varphi 2R2$. Furthermore, the births from $R3$ produce recovered non-pregnant calves, which enter $R1$ giving the additional term $\varphi 3R3$ entering $R1$.

Each of the compartments $S,I,R1,R2$, and $R3$ is subject to a natural death rate $\mu $. However, in the case of $P$, the death rate is increased by $b$ due to the increased mortality of PI animals. Additional births occur in compartments $S,P$, and $R1$. Births from $S$ and $I$ move into the susceptible compartment at rate $\mu $, whereas the births from $P$ occur at a lower rate $\mu \u2212a$ and stay in $P$.

Calculating the change in the total density $N=S+I+R1+R2+R3+P$ results in

Avoiding negative compartment variables $I$, $R2$, $R3$, $P$, which are biologically unfeasible, it follows that

In the proposed model, we aim at keeping the herd size constant. For this purpose, we assume that the reduction in herd density due to $h(I,R2,R3,P)$ is compensated by the introduction of susceptible animals. Furthermore, we assume that the increase in herd density due to $\varphi 3R3$ is compensated by removing animals from the herd regardless of their status. This analysis gives rise to the following set of equations (cf. Fig. 1):

It can be easily seen that the vector field of (3) at the boundary of $(R\u22650)6$ does not point out of $(R\u22650)6$. Therefore, solutions of model (3) are non-negative for all $t\u22650$ if the initial conditions are non-negative. This is an important requirement for an epidemiological model to be meaningful and is not met in the model by Cherry *et al.*^{33}

### B. Equilibria and stability

By defining $F,V$, and $g$ as follows:

we can rewrite model (3) as

where $x=[I,R2,P]$ denotes the disease compartments and $y=[S,R1,R3]$ the disease-free compartments. The model satisfies the following conditions and is, therefore, a well posed epidemiological model according to Chap. 2 of Ref. 34:

$Fi(x,y)\u22650$ for all non-negative $x$ and $y$ and $i=1,2,3$.

$Fi(0,y)=0$ and $Vi(0,y)=0$ for all non-negative $y$ and $i=1,2,3$, i.e., the disease-free set $(0,y)$ is an invariant set.

$Vi(x,y)\u22640$ if $xi=0$ for $i=1,2,3$.

$\u2211i=1nVi(x,y)\u22650$ for all non-negative $x$ and $y$.

The disease-free system $dyidt=gi(0,y)$ has a unique asymptotically stable equilibrium $y0$. Considering assumption (ii), this is an equilibrium of whole system (5) and is called the disease-free equilibrium.

Linearization around the disease-free equilibrium leads to

where $F$ and $V$ are equal to

The dominant eigenvalue of $FV\u22121$ equals the basic reproduction number $R0$ and determines the stability of the disease-free equilibrium.^{35} The disease-free equilibrium is locally asymptotically stable if $R0<1$

For the parameters in Table I, the basic reproduction number equals $7.035$. This result is in agreement with the unstable disease-free equilibrium in Fig. 2.

A successful control measure may be achieved by choosing the removal rate for PI animals above a critical value using the expression for $R0$. A necessary condition to achieve $R0<1$ by increasing the removal rate of PI animals is

The role of TI animals as the source of infection is not entirely clear.^{36} Therefore, in Fig. 3, the dependence of $R0$ on $b$ and $\beta I$ is shown. To gain an insight, we calculate the necessary additional removal rate for PI animals in the case of $\beta I=0$ and find $b>0.098/day$. Thus, we end up with a total removal rate from compartment $P$ of $(\mu +b)=0.099/day$. This means that PI animals should be removed from the herd before reaching an age of $1/(\mu +b)\u224810$ days.

It is also possible to prove the global stability of the disease-free equilibrium for $R0<1$. For this purpose, we consider the compartmental model rewritten as follows:

where $A$ is given by

and $f^$ equals

Then, the following theorem as proven in Chap. 2 of Ref. 34 holds: *If A is a non-singular M-matrix and $f^\u22650$, then the disease-free equilibrium is globally asymptotically stable*. Obviously, $f^\u22650$ is true. An M-matrix can be defined as Z-matrix with eigenvalues whose real parts are non-negative and a Z-matrix is a matrix whose off-diagonal entries are less than or equal to zero. A sufficient condition for a non-singular Z-matrix to be a M-matrix is that it has all non-negative column sums.^{37} This tells us that $V$ is a non-singular M-matrix. Next, we use the following proposition: *if $F$ is non-negative and $V$ is a non-singular M-matrix, then $\rho (FV\u22121)<1$ if and only if all eigenvalues of $(V\u2212F)$ have positive real parts*.^{34} Therefore, we conclude that in the case of $\rho (FV\u22121)<1$ the Z-matrix $A$ is a non-singular M-matrix and the disease-free equilibrium is globally asymptotically stable.

As seen in Fig. 2 and in accordance with the basic reproduction number equal to $7.035$, model (3) approaches the endemic equilibrium. Numerical calculation of the endemic equilibrium results in

The predicted PI animal prevalence of 1.5% is in good agreement with the prevalence of 1.4% found in Danish dairy herds.^{38}

### C. Removal of PI hosts

To simulate a situation where all PI calves are removed from the herd at birth the differential equation for $P$ and the PI transmission coefficient are removed from model (3). In addition, the term $\theta \varphi 2R2$ is removed from the equation for susceptible hosts to adjust for constant total host density. This model assumes that all calves are tested for BVD at birth and that the test is 100% sensitive and specific,

The basic reproduction number of this model equals

An example of the behavior of model (14) is shown in Fig. 4, where the introduction of TI animals in a herd of susceptible animals is simulated. After the initial outbreak, the fraction of recovered hosts decreases due to the removal rate and the fraction of susceptible animals increases. The disease can spread again as soon as the density of susceptible hosts is large enough. The peak prevalence of the subsequent outbreak is reduced due to the presence of recovered animals at the beginning of the outbreak. This behavior results in a damped oscillations of recurrent outbreaks with decreasing peak prevalence, which ultimately approaches an endemic equilibrium. A necessary condition for an outbreak is the presence of TI animals in the herd. This condition is given in Fig. 4 because the density of TI hosts does not reach zero between the outbreaks. However, compartmental models are not a good approximation if the number of hosts in a compartment is low. In a more realistic model, the density of TI animals could reach zero after an outbreak and no subsequent outbreak would be possible. Nevertheless, the simulation of Eq. (14) provides important information about the time required after an outbreak before the herd is susceptible to an outbreak again. This time is mainly determined by $\mu $, which is a measure for the herd turnover rate and can be adjusted by the farmer. Figure 5 visualizes this dependence assuming that the timescale of the oscillating behavior is determined by the complex part of the eigenvalues of the endemic equilibrium if the displacement from the equilibrium is not too large.

## III. STOCHASTIC TRANSMISSION COEFFICIENT

Replacing the transmission coefficient $\beta I$ in the case where all PI calves are removed from the herd at birth with a stochastic transmission coefficient,

where $\xi $ is Gaussian white noise with zero mean and unity variance, leads to the following stochastic compartmental model:

To keep the probability for a negative transmission coefficient negligible the maximum noise intensity $\sigma $ is set to $0.05/day$. At $\sigma =0.05/day$, the probability for negative transmission coefficient equals 0.4%. Figure 6 compares the deterministic BVD model without PI animals with an example path of the BVD model with a stochastic transmission coefficient. In contrast to the deterministic version, the stochastic version is characterized by sustained oscillations. Figure 7 shows the power spectral density for multiple noise intensities and initial conditions equal to the equilibrium value of the deterministic model. The clear peaks in the power spectral densities indicate nearly regular oscillations. The peak positions are equal to the oscillation timescale predicted by the complex eigenvalue of the endemic equilibrium.

To gain an understanding of Fig. 7, we try to derive an expression for the power spectral density based on some simplifications. Since the total density is constant, it is sufficient to analyze the four-dimensional model. Near the deterministic endemic equilibrium equation (17) can be approximated by linearizing the drift coefficient around the endemic equilibrium and replacing $S$ and $I$ in the stochastic term by $S\u22c6=1\u2212I\u22c6\u2212R1\u22c6\u2212R2\u22c6\u2212R3\u22c6$ and $I\u22c6$,

where $x$ and $r$ equal

The endemic equilibrium of the deterministic model without PI hosts calculated numerically equals

The Jacobian matrix evaluated at the endemic equilibrium results in

where $\Phi $ and $\Psi $ equal

Calculating the Fourier transform of Eq. (18) results in

Bringing all terms to the right side leads to

Next, we perform a matrix multiplication from the left side and look at the resulting equation for $I^(f)$,

Finally, we calculate the expected value of the squared modulus of $I^(f)$ to obtain the power spectral density,

It follows that the power spectral density is proportional to the square of the noise intensity. This characteristic is confirmed numerically in the inset plot of Fig. 7, which indicates an approximately constant full width at half maximum of the peaks in the power spectral densities within the investigated noise level. This can be explained by the considered level of noise intensities, which are chosen to keep the model in a biologically plausible range. In other words, increasing the noise intensity within the investigated range increases the amplitude of the oscillation but has no effect on its regularity. This resembles earlier studies on the van der Pol oscillator subject to white noise^{39} and noise-induced oscillators in lasers.^{40}

## IV. DISCUSSION AND OUTLOOK

Modeling the complex spreading dynamics of bovine viral diarrhea remains a challenging task. Based on previous research, we have developed a well posed epidemiological compartmental model that simulates the spreading dynamics within a herd with constant size. The predicted endemic equilibrium of 1.5% is in good agreement with the prevalence of 1.4% found in Danish dairy herds.^{38} The basic reproduction number indicates that increasing the removal rate of PI hosts is a successful control strategy if the transmission coefficient from TI animals is small. This finding is in agreement with the fact that the removal of PI animals is the central component of several effective control strategies.^{41} The removal of PI animals was found to be effective in other simulations as well.^{32} The removal of PI hosts at birth in the deterministic compartmental model results in recurrent outbreaks with decreasing peak prevalence.

To overcome some limitations of the deterministic compartmental model, we have studied a stochastic version that includes randomness in the transmission coefficient. In contrast to the deterministic compartmental model, the model with stochastic transmission coefficient shows sustained oscillations in the case where all PI hosts are removed at birth. Noise-sustained oscillations have been found in many stochastic systems including epidemiological models.^{42,43} In our case, the power spectral density of the sustained oscillations is within the investigated, biologically meaningful noise level proportional to the square of the noise intensity. This is in contrast to the well-known phenomenon of coherence resonance where the coherence of the noise-induced oscillations is maximal for a certain noise intensity.^{44,45} Additional effects might be observed for larger noise intensities such as stochastic bifurcation, which would result in narrower peaks in the power spectrum as known from the van der Pol oscillator.^{46}

Our results suggest many fruitful avenues for future research. The effect of various control strategies could be explored as well as the effect of including vaccination in the compartmental model. Since spatial heterogeneity is highly important in host populations, developing a model involving spatial structure may be of interest. Furthermore, deriving the basic reproduction number in the case of stochastic models could be helpful. Agent-based modeling may be a useful approach to study the transmission dynamics (see, for example, Ref. 30). Agent-based models might underpin explanations of spatial heterogeneity and network interactions in the spreading dynamics of BVD. Furthermore, the deterministic model developed here could be included in more comprehensive models to study the within and between-herd infection dynamics.

## ACKNOWLEDGMENTS

P.H. and K.L. acknowledge support by Deutsche Forschungsgemeinschaft within the framework of Collaborative Research Center 910.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study. All equations and parameters are provided in the text (cf. Table I). The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*et al.*, “

*et al.*, “