Deterministic and stochastic effects in spreading dynamics: A case study of bovine viral diarrhea.

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.


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 Chaos ARTICLE scitation.org/journal/cha persistently viremic for life. 7 Abortions and teratogenic effects can result from infection during approximately the first 150 days. [8][9][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][22][23][24][25] whereas others focus on spreading dynamics. [26][27][28][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][32][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.

A. Model development
The model, schematically shown in Fig. 1, is now described step-by-step. The unit of the compartment variables is hosts/km 2 . 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 p 1 , (2) animals pregnant 1-150 days with fraction p 2 , and (3) animals pregnant 151-280 days with fraction p 3 and p 1 + p 2 + p 3 = 1. • I denotes the fraction of transiently infected (TI) animals.
• P represents the fraction of persistently infected (PI) animals.
• R 1 describes the fraction of recovered, non-pregnant animals.
• R 2 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 R 1 compartment. • R 3 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 R 1 compartment.
Including the above-mentioned fractions p 1 , p 2 , p 3 , 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. 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 (β I I + β P P)S. TI animals move to one of R 1 , R 2 , or R 3 depending on their pregnancy status at the time of infection, i.e., whether they were part of the p 1 S, p 2 S, or p 3 S susceptible subgroups, respectively, giving rise to the corresponding γ p 1 I, γ p 2 I, and γ p 3 I terms.
Recovered animals in R 2 give birth with rate φ 2 and, therefore, move into the non-pregnant recovered compartment R 1 . Likewise, recovered animals in R 3 give birth with rate φ 3 and therefore move into the non-pregnant recovered compartment R 1 . Additionally, births from the R 2 compartment become persistently infected and move into the P compartment; however, due to infection, not all calves survive and the number of births φ 2 R 2 is reduced by the factor θ , giving the term θ φ 2 R 2 . Furthermore, the births from R 3 produce recovered non-pregnant calves, which enter R 1 giving the additional term φ 3 R 3 entering R 1 .
Each of the compartments S, I, R 1 , R 2 , and R 3 is subject to a natural death rate µ. 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 R 1 . Births from S and I move into the susceptible compartment at rate µ, whereas the births from P occur at a lower rate µ − a and stay in P.
Calculating the change in the total density N = S + I + R 1 Avoiding negative compartment variables I, R 2 , R 3 , 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, R 2 , R 3 , P) is compensated by the introduction of susceptible animals. Furthermore, we assume that the increase in herd density due to φ 3 R 3 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 ≥0 ) 6 does not point out of (R ≥0 ) 6 . Therefore, solutions of model (3) are non-negative for all t ≥ 0 if the initial conditions are nonnegative. This is an important requirement for an epidemiological model to be meaningful and is not met in the model by Cherry et al. 33 Figure 2 shows the behavior of model (3) after the introduction of one PI animal per km 2 into a herd of susceptible animals with a host density of 67 animals/km 2 . The outbreak of disease is followed by an approach to the endemic equilibrium.

B. Equilibria and stability
By defining F , V, and g as follows: we can rewrite model (3) as (v) The disease-free system dy i dt = g i (0, y) has a unique asymptotically stable equilibrium y 0 . 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 −1 equals the basic reproduction number R 0 and determines the stability of the disease-free equilibrium. 35 The disease-free equilibrium is locally asymptotically stable if R 0 < 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 R 0 . A necessary condition to achieve R 0 < 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 R 0 on b and β I is shown. To gain an insight, we calculate the necessary additional removal rate for PI animals in the case of β I = 0 and find b > 0.098/day. Thus, we end up with a total removal rate from compartment P of (µ + b) = 0.099/day. This means that PI animals should be removed from the herd before reaching an age of 1/(µ + b) ≈ 10 days. It is also possible to prove the global stability of the diseasefree equilibrium for R 0 < 1. For this purpose, we consider the compartmental model rewritten as follows:  Table I.

ARTICLE scitation.org/journal/cha
where A is given by andf equalŝ Then, the following theorem as proven in Chap. 2

of Ref. 34 holds:
If A is a non-singular M-matrix andf ≥ 0, then the disease-free equilibrium is globally asymptotically stable. Obviously,f ≥ 0 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 offdiagonal 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 nonnegative and V is a non-singular M-matrix, then ρ(FV −1 ) < 1 if and only if all eigenvalues of (V − F) have positive real parts. 34 Therefore, we conclude that in the case of ρ(FV −1 ) < 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)  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 θ φ 2 R 2 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 µ, 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.  Table I. The calculation of the period is based on the complex part of the eigenvalues of the endemic equilibrium.

III. STOCHASTIC TRANSMISSION COEFFICIENT
Replacing the transmission coefficient β I in the case where all PI calves are removed from the herd at birth with a stochastic transmission coefficient, where ξ 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 σ is set to 0.05/day. At σ = 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 fourdimensional 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 where x and r equal x = I, R 1 , R 2 , R 3 , r = σ I S , 0, 0, 0 .

FIG. 7.
Power spectral densities of model (17) for TI hosts at different noise intensities σ (color bar) calculated from 300 time series with a simulation period of 10 4 days and initial conditions at the deterministic endemic equilibrium. The dashed gray line represents the oscillation frequency predicted by the deterministic model according to Fig. 5. Inset: dependence of the full width at half maximum of the peaks (FWHM in 10 −4 /day) on the noise intensity.

ARTICLE scitation.org/journal/cha
The endemic equilibrium of the deterministic model without PI hosts calculated numerically equals S = 43.37%, I =0.89%, R 1 =52.70%, The Jacobian matrix evaluated at the endemic equilibrium results in where and equal Calculating the Fourier transform of Eq.
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Î(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.

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.