Stochastic Modelling and Analysis of Within-Farm Highly Pathogenic Avian Influenza Dynamics in Dairy Cattle

1. Introduction

Avian influenza (AI), commonly referred to as bird flu, is a highly contagious viral disease caused by influenza A viruses belonging to the family Orthomyxoviridae and the genus Influenza A [1,2]. These viruses primarily infect avian species, especially wild aquatic birds, which act as the natural reservoir and contribute to global transmission through migratory movements [3,4]. While wild birds are the primary host, spillover into domestic poultry can lead to widespread outbreaks, causing severe economic and ecological consequences. Avian influenza viruses are classified into two categories: low pathogenic avian influenza (LPAI), which causes mild or subclinical infections, and highly pathogenic avian influenza (HPAI), which results in severe systemic disease and high mortality rates [5].

Globally, HPAI H5N1 subtype remains one of the most concerning strains due to its ability to infect multiple species and its persistence in the environment. Currently, only New Zealand, Australia, and Antarctica remain free from HPAI H5N1 [6,7], outbreaks in other parts of the world highlight the virus’s capacity for rapid adaptation and cross-species transmission. This cross-species potential has led to heightened global vigilance in monitoring wild bird populations, their environment and domestic animal interfaces.

The recent emergence of HPAI H5N1 in non-avian hosts, particularly in dairy cattle, marks a significant epidemiological shift [8]. In 2024, the United States confirmed several outbreaks of H5N1 in dairy herds, with detections in nasal swabs and unpasteurized milk, suggesting viral replication within the mammary gland and possible environmental contamination through the milking system [5,9,10,11,12]. The initial cases were detected in Texas and subsequently in Kansas, Michigan, and New Mexico. These findings challenge the long-held assumption that cattle are resistant to avian influenza infection and highlight new pathways for viral maintenance outside avian populations [13,14]. The transmission was linked to direct cow-to-cow contact, contaminated milking equipment, human-assisted spread, and potential exposure to infected wild birds or poultry. Clinically, affected cattle exhibited reduced milk production, fever, nasal discharge, and, in many cases, mild or asymptomatic presentations, complicating detection [15,16].

The incubation phase of HPAI in cattle can last up to two weeks; during this window, the movement of asymptomatic carriers between farms creates a risk of further geographic dissemination [17]. Although H5N1 is highly lethal in avian species, mortality and culling rates in cattle have historically remained low, typically staying below 2% [18]. While most infected cattle recover with supportive treatment, the virus continues to present substantial economic and biosecurity challenges [19]. Immunity gained from prior infection is expected to be partial or transient—particularly given viral evolution—making lifelong protection improbable [20]. Consequently, without rigorous surveillance and biosecurity protocols, recurrent outbreaks remain a distinct possibility [21].

Avian influenza transmission occurs through various pathways, including direct interaction with infected hosts, inhalation or ingestion of viral particulates, contact with contaminated fomites, and exposure to environmental reservoirs [21,22]. While bird-to-bird transmission is well documented, the identification of HPAI in mammals raises concerns regarding environmental persistence and potential adaptation to mammalian physiology. Occupational groups such as veterinarians, poultry and dairy workers, and animal handlers face the highest risk of exposure. In birds, clinical signs often include neurological deficits—such as tremors, torticollis, and ataxia—alongside respiratory distress and diarrhoea, with severe cases frequently resulting in sudden death. Although human infections are rare, they can manifest as fever, muscle pain, and respiratory issues, potentially progressing to pneumonia or multi-organ failure [16].

The emergence of HPAI H5N1 in cattle represents a major shift in the virus’s epidemiological behavior, with significant consequences for animal health, public health, and the agricultural industry [5,9]. Effective control strategies require a detailed understanding of how the virus spreads within herds–especially the relative contribution of direct animal-to-animal transmission compared with environmental pathways, and the conditions that enable outbreaks to continue [21,23,24]. Mathematical modelling offers a powerful framework for investigating these issues by quantifying distinct transmission modes, pinpointing key epidemiological thresholds, and assessing the likely impact of different control options [25,26,27]. As the epidemiology of HPAI evolves, there is a pressing need for mathematical models that explicitly incorporate both direct and environmental transmission processes. These models are essential for estimating infection burdens, determining critical control levers, and comparing the performance of intervention strategies across a range of epidemiological settings. Recently, the transmission dynamics of Avian influenza in cattle population in the United States has been studied using mathematical models in [14].

Motivated by the recent HPAI H5N1 outbreak in dairy cattle in the US [12,14], we develop and analyse a stochastic within-farm model for HPAI transmission. The model incorporates both symptomatic and asymptomatic cattle, as well as virus contamination in the environment. Using standard stochastic differential equation theory, we establish the positivity and boundedness of solutions and derive conditions for the extinction and persistence of HPAI. Numerical experiments are carried out to support and illustrate the theoretical findings.

The remainder of the paper is structured as follows: Section 2 presents the formulation of the stochastic HPAI model. Section 3 is devoted to the qualitative analysis of the model. Numerical simulations are reported in Section 4. Finally, the conclusion and study limitations are given in Section 5 and Section 6, respectively.

2. Model Formulation

The following key assumptions guide the model formulation:

• The cattle herd is treated as an open population: animals enter through recruitment, purchase, or births at a constant rate $\Lambda$, and leave through natural mortality at rate $\mu$ or disease-induced mortality at rate d. This setup reflects realistic herd management practices rather than a fully closed herd.
• Two infection pathways are considered: (a) direct transmission, when susceptible cattle have contact with infectious animals at the rate ${\beta }_a$ and ${\beta }_s$, and (b) indirect transmission at rate ${\beta }_B$, when cattle encounter virus particles present in the environment (such as on contaminated surfaces, bedding, or milking equipment).
• Infected animals may develop clinical signs or remain asymptomatic. Both groups shed virus, but asymptomatic cattle typically go unnoticed, prolonging their infectious period and helping to maintain transmission within the herd.
• Once infected animals are identified and removed from the herd through culling, separation, slaughter, or other means, they are assumed no longer to contribute to transmission and are placed in the removed class (R).
• The environmental compartment (B) represents viral particles deposited into the environment by infected cattle. This environmental burden declines at the constant rate $\epsilon$, reflecting loss of infectivity or degradation under natural conditions (e.g., temperature, sunlight, disinfectants).
• Infection arising from environmental exposure is modelled using a saturating function, so that the risk of infection increases with contamination but levels off once exposure becomes very high.

Based on the assumptions, a mathematical model named $SE{I}_s{I}_{a}R-B$ for HPAI was developed where total cattle population is divided into five compartments: Susceptible (S), Exposed (E), Symptomatic ($I_s$), Asymptomatic ($I_a$), Removed (R), and contaminated environment (B). The model equations are written using a system of ordinary differential equations (ODEs) and given by

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -{\beta }_{s}S{\displaystyle \frac{I_s}{N}}-{\beta }_{a}S{\displaystyle \frac{I_a}{N}}-{\beta }_{B}S{\displaystyle \frac{B}{K+B}}-\mu S,\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& ={\beta }_{s}S{\displaystyle \frac{I_s}{N}}+{\beta }_{a}S{\displaystyle \frac{I_a}{N}}+{\beta }_{B}S{\displaystyle \frac{B}{K+B}}-\sigma E-\mu E,\hfill \\ \hfill {\displaystyle \frac{d{I}_s}{dt}}& =\nu \sigma E-(\mu +d+\gamma )I_s,\hfill \\ \hfill {\displaystyle \frac{d{I}_a}{dt}}& =(1-\nu )\sigma E-(\mu +d+\delta )I_a,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =\gamma I_s+\delta I_a-\mu R,\hfill \\ \hfill {\displaystyle \frac{dB}{dt}}& ={\omega }_s{I}_s+{\omega }_a{I}_a-\epsilon B.\hfill \end{array}$$

(1)

with initial conditions

$$S\left(0\right)\ge 0,E\left(0\right)\ge 0,I_s\left(0\right)\ge 0,I_a\left(0\right)\ge 0,R\left(0\right)\ge 0,B\left(0\right)\ge 0.$$

The total population at any time t satisfies

$$N\left(t\right)=S+E+I_s+I_a+R.$$

A schematic diagram of the model structure is presented in Figure 1, illustrating the flow of individuals between compartments and the interactions between cattle and their environment. Each arrow denotes a transition governed by a specific rate parameter, such as infection, progression, recovery, or decay.

To capture the inherent randomness and environmental variability influencing the transmission dynamics of avian influenza (AI), we extend the deterministic compartmental framework to a stochastic differential equation (SDE) system by incorporating multiplicative white noise perturbations in each epidemiological class except recovery class as it will not affect the dynamics of infected class [28]. The deterministic model assumes perfect homogeneity in disease transmission and parameter constancy over time; however, in natural ecological settings such as cattle migration intensity, pathogen survival in the environment, and host immune variability can significantly alter infection rates and removal processes [29,30]. The stochastic version of the model (2) therefore provides a more realistic representation of AI dynamics by allowing such uncertainties to enter the system as diffusion terms proportional to the state variables.

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\Lambda -{\beta }_{s}S{\displaystyle \frac{I_s}{N}}-{\beta }_{a}S{\displaystyle \frac{I_a}{N}}-{\beta }_{B}S{\displaystyle \frac{B}{K+B}}-\mu S+{\sigma }_{S}S{\displaystyle \frac{d{W}_S}{dt}}\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& ={\beta }_{s}S{\displaystyle \frac{I_s}{N}}+{\beta }_{a}S{\displaystyle \frac{I_a}{N}}+{\beta }_{B}S{\displaystyle \frac{B}{K+B}}-\sigma E-\mu E+{\sigma }_{e}E{\displaystyle \frac{d{W}_e}{dt}}\hfill \\ \hfill {\displaystyle \frac{d{I}_s}{dt}}& =\nu \sigma E-(\mu +d+\gamma )I_s+{\sigma }_s{I}_s{\displaystyle \frac{d{W}_s}{dt}}\hfill \\ \hfill {\displaystyle \frac{d{I}_a}{dt}}& =(1-\nu )\sigma E-(\mu +d+\delta )I_a+{\sigma }_a{I}_a{\displaystyle \frac{d{W}_a}{dt}}\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =\gamma I_s+\delta I_a-\mu R\hfill \\ \hfill {\displaystyle \frac{dB}{dt}}& ={\omega }_s{I}_s+{\omega }_a{I}_a-\epsilon B+{\sigma }_{B}B{\displaystyle \frac{d{W}_B}{dt}}\hfill \end{array}$$

(2)

Each compartment is influenced by an independent Wiener process $W_i\left(t\right),i=S,e,s,a,B$ representing the cumulative effect of random perturbations acting on that class. The term ${\sigma }_{S}Sd{W}_S\left(t\right)$ captures random fluctuations in the recruitment and natural death processes of susceptible hosts due to unpredictable environmental events that changes host availability or contact rates. The exposed class includes the diffusion term ${\sigma }_{e}Ed{W}_e\left(t\right)$, which reflects uncertainty in the latent period, for instance, due to temperature-dependent virus incubation to infected excreta. The stochastic perturbation ${\sigma }_s{I}_{s}d{W}_s\left(t\right)$ in the symptomatic infected class models random variability in transmission and recovery among visibly infected individuals, possibly caused by variations in immune response. Similarly, ${\sigma }_a{I}_{a}d{W}_a\left(t\right)$ in the asymptomatic infected compartment accounts for random changes in subclinical infection dynamics, such as variations in viral shedding or host contact behavior that are difficult to monitor in the field. The environmental virus reservoir $B\left(t\right)$ is perturbed by ${\sigma }_{B}Bd{W}_B\left(t\right)$, representing unpredictable changes in viral persistence and decay in the environment, affected by climatic and ecological factors [31].

Each noise term is multiplicative which shows that the magnitude of random fluctuations scales with the population level of the corresponding compartment. This is biologically feasible because the stochastic influence vanishes when the compartment population is zero and increases with population density. Such a formulation guarantees positivity of solutions and preserves the nonnegativity of all compartments [32]. The stochastic framework thus enables investigation of how random environmental shocks and intrinsic variability can drive disease persistence, extinction, or oscillations even when deterministic conditions predict eradication. The parameters used in the model and their biological interpretation are summarised in Table 1.

3. Analysis of Model

3.1. Positivity and Boundedness

Consider the stochastic AI system in Itô form

$$\left.\begin{array}{cc}\hfill dS& =\left[\Lambda -{\beta }_s\frac{S{I}_s}{N}-{\beta }_a\frac{S{I}_a}{N}-{\beta }_B\frac{SB}{K+B}-\mu S\right]dt+{\sigma }_{S}Sd{W}_S\left(t\right)\hfill \\ \hfill dE& =\left[{\beta }_s\frac{S{I}_s}{N}+{\beta }_a\frac{S{I}_a}{N}+{\beta }_B\frac{SB}{K+B}-\sigma E-\mu E\right]dt+{\sigma }_{E}Ed{W}_e\left(t\right)\hfill \\ \hfill d{I}_s& =\left[\nu \sigma E-(\mu +d+\gamma )I_s\right]dt+{\sigma }_s{I}_{s}d{W}_s\left(t\right)\hfill \\ \hfill d{I}_a& =\left[(1-\nu )\sigma E-(\delta +\mu +d)I_a\right]dt+{\sigma }_a{I}_{a}d{W}_a\left(t\right)\hfill \\ \hfill dR& =(\gamma I_s+\delta I_a-\mu R)dt\hfill \\ \hfill dB& =({\omega }_s{I}_s+{\omega }_a{I}_a-\epsilon B)dt+{\sigma }_{B}Bd{W}_B\left(t\right)\hfill \end{array}\right\}$$

(3)

where $W_i\left(t\right)$ are independent standard Brownian motions.

3.1.1. Positivity

Assuming the standard hypotheses that all parameters are nonnegative and the initial state

$$(S_0,E_0,I_{s,0},I_{a,0},R_0,B_0)\in {\mathbb{R}}_{+}^6$$

is nonnegative a.s. Then the Itô SDE system admits a unique global strong solution $X_t=(S_t,E_t,I_{s,t},I_{a,t},R_t,B_t)$ with

$$X_t\in {\mathbb{R}}_{+}^6\mathrm{for}\mathrm{all}t\ge 0\mathrm{a}.\mathrm{s}.,$$

and the first moments satisfy uniform-in-time bounds and hence the solution is bounded.

Proof. 

Since the diffusion terms are proportional to the state variables (${\sigma }_i{X}_{i}d{W}_i$), we can apply Itô’s formula to each class of model (3) for every strictly positive component $X_t$. Writing SDE in compact form as

$$d{X}_t=f_X\left(X_t\right)dt+{\sigma }_X{X}_{t}d{W}_t$$

(4)

where $f_X(·)$ denotes the drift of that component and is locally Lipschitz. Apply Itô’s formula to $ln{X}_t$ while $X_t>0$:

$$d(ln{X}_t)=\left({\displaystyle \frac{f_X\left(X_t\right)}{X_t}}-\frac{1}{2}{\sigma }_X^2\right)dt+{\sigma }_{X}d{W}_t.$$

(5)

The right-hand side of (5) is finite for all finite times before explosion because $f_X$ is locally bounded on compacts and $X_t>0$ on that interval.

This means, for any bounded subset $P\subset {\mathbb{R}}_{+}^6$, there exists a constant $C_P>0$ such that

$$|f_X\left(x\right)|\le C_P,\forall x\in P.$$

Thus, whenever the state variables $(S_t,E_t,I_{s,t},I_{a,t},R_t,B_t)$ remain within a finite range, the drift term $f_X$ cannot become infinite. Since all drift terms in the stochastic AI model are composed of polynomial and rational expressions with positive denominators, they are locally bounded on the positive orthant ${\mathbb{R}}_{+}^6$. This ensures that the stochastic integrals ${\int }_0^t{f}_X\left(X_z\right)dz$ are well defined and finite for all finite t.

Now, When applying Itô’s formula to $ln{X}_t$, it is necessary that $X_t>0$, so that $ln\left(X_t\right)$ is well defined. Therefore, the analysis is carried out on a maximal time interval $[0,\tau )$ for which $X_t>0$ almost surely. If we assume, for contradiction, that $X_t$ reaches zero at a finite time $\tau$, then $ln\left(X_t\right)\to -\infty$ as $t\to \tau$. However, since $f_X$ is locally bounded and ${\sigma }_X$ is finite, the right–hand side of the Itô differential (5) remains finite on $[0,\tau )$, leading to a contradiction. Hence $X_t$ cannot reach zero in finite time and remains positive almost surely.

Also, the recovered class $R_t$ has no diffusion term and its drift is $\gamma I_{s,t}+\delta I_{a,t}-\mu R_t$. Since $I_{s,t},I_{a,t}\ge 0$ a.s., standard ODE comparison shows $R_t$ cannot become negative even if $R_0\ge 0$.

Hence all components remain strictly positive whenever their initial value is positive for all $t\ge 0$ almost surely. □

3.1.2. Boundedness

We have

$$N_t:=S_t+E_t+I_{s,t}+I_{a,t}+R_t.$$

Summing the SDEs (3) for $S_t,E_t,I_{s,t},I_{a,t},R_t$ yields

$$d{N}_t=\left(\Lambda -\mu N_t-d{I}_{s,t}\right)dt+{\sigma }_S{S}_{t}d{W}_S+{\sigma }_e{E}_{t}d{W}_e+{\sigma }_s{I}_{s,t}d{W}_s+{\sigma }_a{I}_{a,t}d{W}_a.$$

Taking expectation and using $\mathbb{E}\left[d{W}_i\right]=0$, we get

$${\displaystyle \frac{d}{dt}}\mathbb{E}\left[N_t\right]=\Lambda -\mu \mathbb{E}\left[N_t\right]-d\mathbb{E}\left[I_{s,t}\right]\le \Lambda -\mu \mathbb{E}\left[N_t\right]$$

Solving this linear differential inequality, we get

$$\mathbb{E}\left[N_t\right]\le \mathbb{E}\left[N_0\right]e^{-\mu t}+{\displaystyle \frac{\Lambda }{\mu }}\left(1-e^{-\mu t}\right)\le max\{\mathbb{E}\left[N_0\right],\Lambda /\mu \},$$

so the host total is uniformly bounded in expectation.

Now, we consider the environmental compartment $B_t$. Using SDE from (3),

$${\displaystyle \frac{d}{dt}}\mathbb{E}\left[B_t\right]=\mathbb{E}[{\omega }_s{I}_{s,t}+{\omega }_a{I}_{a,t}-\epsilon B_t]$$

Using the bound $I_{s,t}+I_{a,t}\le N_t$ and let ${\omega }_{max}=max\{{\omega }_s,{\omega }_a\}$,

$${\displaystyle \frac{d}{dt}}\mathbb{E}\left[B_t\right]\le {\omega }_{max}\mathbb{E}\left[N_t\right]-\epsilon \mathbb{E}\left[B_t\right]$$

Since $\mathbb{E}\left[N_t\right]$ is uniformly bounded, this is a linear inhomogeneous differential inequality for $\mathbb{E}\left[B_t\right]$. Solving using variation of constants gives

$$\mathbb{E}\left[B_t\right]\le \mathbb{E}\left[B_0\right]e^{-\epsilon t}+{\displaystyle \frac{{\omega }_{max}{sup}_{s\ge 0}\mathbb{E}\left[N_s\right]}{\epsilon }}\left(1-e^{-\epsilon t}\right),$$

hence ${sup}_{t\ge 0}\mathbb{E}\left[B_t\right]<\infty$

Thus, every component of model (3) has a finite first moment uniformly in time. Using standard SDE theory [33], bounded first moments prevent finite-time explosion and thus the local solution extends to all $t\ge 0$.

3.2. Equilibrium Analysis

3.2.1. Disease-Free Equilibrium (DFE)

At the disease-free equilibrium (DFE), the infection is completely absent from the cattle population, meaning that no individuals exist in the exposed or infectious compartments. To obtain this equilibrium, we set all time derivatives in the system to zero and assume that no disease transmission occurs ($E=I_s=I_a=R=B=0$). Under these conditions, the number of susceptible cattle remains constant, determined only by the balance between recruitment and natural death. The noise term in stochastc model does not contribute to equilibrium itself because its expectation is zero ($\mathbb{E}(d{W}_i\left(t\right)=0$). However, noise does affect the stability and long-term behaviour around that equilibrium [31]. At the disease–free equilibrium all infected compartments vanish,

so set $E^{*}=I_s^{*}=I_a^{*}=R^{*}=B^{*}=0$. Thus we have

$$0=\Lambda -\mu S^{*}\Rightarrow S^{*}={\displaystyle \frac{\Lambda }{\mu }}.$$

Hence the DFE is

$$\mathrm{DFE}=\left(S^{*},E^{*},I_s^{*},I_a^{*},R^{*},B^{*}\right)=\left(\frac{\Lambda }{\mu },0,0,0,0,0\right).$$

3.2.2. Basic Reproduction Number ${\mathcal{R}}_0$

The basic reproduction number, denoted as ${\mathcal{R}}_0$, represents the average number of secondary infections produced by a single infectious cattle introduced into a completely susceptible population [26]. It serves as a key threshold parameter that determines whether an infection can invade and persist within a population.

We derive ${\mathcal{R}}_0$ using the next–generation matrix denoted by ${\psi }_{NGM}$. At the DFE we have $N^{*}=S^{*}$, hence ${\displaystyle \frac{S^{*}}{N^{*}}}=1$. Therefore the linearised new–infection matrix F at DFE is

$$F=\left(\begin{array}{cccc}0& {\beta }_s& {\beta }_a& {\displaystyle \frac{{\beta }_B{S}^{*}}{K}}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$

The transition matrix V describes the flow of individuals out of infected compartments due to recovery, removal or death. The transition matrix V is given by

$$V=\left(\begin{array}{cccc}\sigma +\mu & 0& 0& 0\\ -\nu \sigma & \mu +d+\gamma & 0& 0\\ (\nu -1)\sigma & 0& \mu +\delta +d& 0\\ 0& -{\omega }_s& -{\omega }_a& \epsilon \end{array}\right).$$

Rows correspond to the ordering $(E,I_s,I_a,B)$. By construction the diagonal entries are the rates leaving each infected class.

The next generation matrix ${\psi }_{NGM}=F{V}^{-1}$ and the basic reproduction number is the spectral radius of the next–generation matrix. Thus,

$${\mathcal{R}}_0=\rho \left({\psi }_{NGM}\right)=\rho \left(F{V}^{-1}\right).$$

We have

$$V^{-1}=\left(\begin{array}{cccc}{\displaystyle \frac{1}{\sigma +\mu }}& 0& 0& 0\\ {\displaystyle \frac{\nu \sigma }{(\sigma +\mu )(\mu +d+\gamma )}}& {\displaystyle \frac{1}{\mu +d+\gamma }}& 0& 0\\ {\displaystyle \frac{(1-\nu )\sigma }{(\sigma +\mu )(\mu +\delta +d)}}& 0& {\displaystyle \frac{1}{\mu +\delta +d}}& 0\\ {\displaystyle \frac{\nu \sigma {\omega }_s}{\epsilon (\sigma +\mu )(\mu +d+\gamma )}}+{\displaystyle \frac{(1-\nu )\sigma {\omega }_a}{\epsilon (\sigma +\mu )(\mu +\delta +d)}}& {\displaystyle \frac{{\omega }_s}{\epsilon (\mu +d+\gamma )}}& {\displaystyle \frac{{\omega }_a}{\epsilon (\mu +\delta +d)}}& {\displaystyle \frac{1}{\epsilon }}\end{array}\right).$$

Thus, we get the basic reproduction number

$${\mathcal{R}}_0={\displaystyle \frac{\sigma }{\sigma +\mu }}\left[{\displaystyle \frac{\nu {\beta }_s}{\mu +d+\gamma }}+{\displaystyle \frac{(1-\nu ){\beta }_a}{\mu +\delta +d}}+{\displaystyle \frac{{\beta }_B}{K\epsilon }}\left({\displaystyle \frac{\nu {\omega }_s}{\mu +d+\gamma }}+{\displaystyle \frac{(1-\nu ){\omega }_a}{\mu +\delta +d}}\right)\right].$$

(6)

The prefactor $\sigma /(\sigma +\mu )$ accounts for the probability that an exposed host survives the latent period and becomes infectious. The first two terms in the bracket are the sum of direct (host-to-host) contributions from symptomatic and asymptomatic infectives. The last two terms in the bracket quantifies the indirect (environmental) transmission that an infected host sheds virus into the environment at rates ${\omega }_s,{\omega }_a$), this virus decays at rate $\epsilon$, and environmental virus produces new infections at rate proportional to ${\beta }_B{S}^{*}/\left(K\epsilon \right)$. The stochastic diffusion terms in the model do not change the location of the DFE (they vanish when $E=I_s=I_a=B=0$), but they do affect stability and sample-path behaviour near the DFE. The stochastic stability analysis near the equilibrium is discussed in next section.

3.2.3. Endemic Equilibrium (EE)

We find the endemic equilibrium of the deterministic part of the system. Set all stochastic terms and time derivatives to zero. When we analyse endemic equilibrium in epidemic models, we mean a steady state of the system where the number of individuals in each compartment remains constant over time. But, Brownian motion has no steady state: its increments are unbounded, so the noise term cannot vanish at a fixed positive level. The equilibrium $(S^{**},E^{**},I_s^{**},I_a^{**},R^{**},B^{**})$ satisfies the algebraic system

$$\begin{array}{cc}\hfill 0& =\Lambda -{\beta }_a{\displaystyle \frac{S^{**}{I}_a^{**}}{N^{**}}}-{\beta }_s{\displaystyle \frac{S^{**}{I}_s^{**}}{N^{**}}}-{\beta }_B{\displaystyle \frac{S^{**}{B}^{**}}{K+B^{**}}}-\mu S^{**},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill 0& ={\beta }_a{\displaystyle \frac{S^{**}{I}_a^{**}}{N^{**}}}+{\beta }_s{\displaystyle \frac{S^{**}{I}_s^{**}}{N^{**}}}+{\beta }_B{\displaystyle \frac{S^{**}{B}^{**}}{K+B^{**}}}-(\sigma +\mu )E^{**},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill 0& =\nu \sigma E^{**}-(\mu +d+\gamma )I_s^{**},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill 0& =(1-\nu )\sigma E^{**}-(\delta +\mu +d)I_a^{**},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill 0& =\gamma I_s^{**}+\delta I_a^{**}-\mu R^{**},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill 0& ={\omega }_s{I}_s^{**}+{\omega }_a{I}_a^{**}-\epsilon B^{**}.\hfill \end{array}$$

(12)

Here $N^{**}=S^{**}+E^{**}+I_s^{**}+I_a^{**}+R^{**}$. We are supposed to find solution with $E^{**}>0$ (and hence $I_s^{**},I_a^{**},B^{**},R^{**}\ge 0$).

First solve the linear equations (9)–(12) in terms of $E^{**}$. From equation (9) and equation (10),

$$\begin{array}{cc}\hfill I_s^{**}& ={\alpha }_s{E}^{**},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill I_a^{**}& ={\alpha }_a{E}^{**}.\hfill \end{array}$$

(14)

where

$${\alpha }_s:={\displaystyle \frac{\nu \sigma }{\mu +d+\gamma }},{\alpha }_a:={\displaystyle \frac{(1-\nu )\sigma }{\delta +\mu +d}}.$$

are the two positive proportionality constants.

From equation (11) we obtain

$$R^{**}={\displaystyle \frac{\gamma I_s^{**}+\delta I_a^{**}}{\mu }}={\displaystyle \frac{\gamma {\alpha }_s+\delta {\alpha }_a}{\mu }}{E}^{**}.$$

Define for convenience

$$\rho :={\displaystyle \frac{\gamma {\alpha }_s+\delta {\alpha }_a}{\mu }},$$

so $R^{**}=\rho E^{**}$. From equation (12) we obtain

$$B^{**}={\displaystyle \frac{{\omega }_s{I}_s^{**}+{\omega }_a{I}_a^{**}}{\epsilon }}={\displaystyle \frac{{\omega }_s{\alpha }_s+{\omega }_a{\alpha }_a}{\epsilon }}{E}^{**}.$$

Define

$$\zeta :={\displaystyle \frac{{\omega }_s{\alpha }_s+{\omega }_a{\alpha }_a}{\epsilon }},$$

so $B^{**}=\zeta E^{**}$. Now using the relations above we can express $N^{**}$ in terms of $S^{**}$ and $E^{**}$.

$$N^{**}=S^{**}+E^{**}+{\alpha }_s{E}^{**}+{\alpha }_a{E}^{**}+\rho E^{**}=S^{**}+c{E}^{**}.$$

where

$$c:=1+{\alpha }_s+{\alpha }_a+\rho =1+{\alpha }_s+{\alpha }_a+{\displaystyle \frac{\gamma {\alpha }_s+\delta {\alpha }_a}{\mu }}.$$

Let us introduce two constants

$$A_1:={\beta }_s{\alpha }_s+{\beta }_a{\alpha }_a,A_2:={\beta }_B\zeta .$$

Here $A_1$ collects the contributions of cattle-to-cattle transmission via symptomatic and asymptomatic infectives, and $A_2$ collects the environmental contribution.

Now use the equation (8) to relate $S^{**}$ and $E^{**}$. Plugging $I_s^{**}={\alpha }_s{E}^{**},I_a^{**}={\alpha }_a{E}^{**},B^{**}=\zeta E^{**}$ into (8) gives

$${\beta }_a{\displaystyle \frac{S^{**}{\alpha }_a{E}^{**}}{N^{**}}}+{\beta }_s{\displaystyle \frac{S^{**}{\alpha }_s{E}^{**}}{N^{**}}}+{\beta }_B{\displaystyle \frac{S^{**}\zeta E^{**}}{K+\zeta E^{**}}}=(\sigma +\mu )E^{**}.$$

Factor $E^{**}$ (and cancel $E^{**}>0$):

$$S^{**}\left({\displaystyle \frac{A_1}{N^{**}}}+{\displaystyle \frac{A_2}{K+\zeta E^{**}}}\right)=\sigma +\mu .$$

Equivalently,

$${\displaystyle \frac{\sigma +\mu }{S^{**}}}={\displaystyle \frac{A_1}{N^{**}}}+{\displaystyle \frac{A_2}{K+\zeta E^{**}}}.$$

This is one relation between $S^{*}$ and $E^{*}$ (recall $N^{*}=S^{*}+c{E}^{*}$). Now, use the equation (7) to express $S^{**}$ by the infection pressure. Let us define the instantaneous infection pressure (per susceptible)

$${\lambda }^{**}:={\displaystyle \frac{{\beta }_s{I}_s^{**}}{N^{**}}}+{\displaystyle \frac{{\beta }_a{I}_a^{**}}{N^{**}}}+{\displaystyle \frac{{\beta }_B{B}^{**}}{K+B^{**}}}.$$

Then equation (7) is

$$0=\Lambda -S^{**}(\mu +{\lambda }^{**}),$$

hence

$$S^{**}={\displaystyle \frac{\Lambda }{\mu +{\lambda }^{**}}}.$$

From equation (8) we have ${\lambda }^{**}{S}^{**}=(\sigma +\mu )E^{**}$ (because the LHS of equation (8) equals ${\lambda }^{**}{S}^{**}$); therefore

$${\lambda }^{**}={\displaystyle \frac{(\sigma +\mu )E^{**}}{S^{**}}}.$$

Substituting equation (16) into equation (17) we eliminate $S^{**}$ and obtain an explicit expression for ${\lambda }^{**}$ as a function of $E^{**}$:

$${\lambda }^{**}={\displaystyle \frac{(\sigma +\mu )E^{**}}{\Lambda /(\mu +{\lambda }^{**})}}⟹{\lambda }^{**}(\Lambda -(\sigma +\mu )E^{**})=(\sigma +\mu )\mu E^{**}.$$

Hence, provided $\Lambda -(\sigma +\mu )E^{**}\ne 0$,

$${\lambda }^{**}={\displaystyle \frac{(\sigma +\mu )\mu E^{**}}{\Lambda -(\sigma +\mu )E^{**}}}.$$

Above equation (18) cancels the trivial solution $E^{**}=0$; the endemic branch requires $E^{**}>0$ and the denominator positive, so a necessary admissibility restriction is

$0<E^{**}<{\displaystyle \frac{\Lambda }{\sigma +\mu }}$.

Now obtain a single scalar equation for $E^{**}$.

We have two expressions for the infection pressure ${\lambda }^{**}$:

(i) By definition and using the proportionalities,

$${\lambda }^{**}=E^{**}\left({\displaystyle \frac{A_1}{N^{**}}}+{\displaystyle \frac{A_2}{K+\zeta E^{**}}}\right),$$

because $I_s^{**}={\alpha }_s{E}^{**},I_a^{**}={\alpha }_a{E}^{**},B^{**}=\zeta E^{**}$ and $N^{**}=S^{**}+c{E}^{**}$.

(ii) By elimination we have equation (18) which gives ${\lambda }^{**}$ solely as a function of $E^{**}$. Equate the equations (18) and (19); cancel $E^{**}$ (since we search for $E^{**}>0$):

$${\displaystyle \frac{(\sigma +\mu )\mu }{\Lambda -(\sigma +\mu )E^{**}}}={\displaystyle \frac{A_1}{N^{**}}}+{\displaystyle \frac{A_2}{K+\zeta E^{**}}}.$$

We have $N^{**}=S^{**}+c{E}^{**}$ and $S^{**}=\Lambda /(\mu +{\lambda }^{**})$ with ${\lambda }^{**}$ given by equation (18). Thus the left-hand side is an explicit function of $E^{**}$, and the right-hand side is also an explicit function of $E^{**}$ after substitution. Equation (20) is the required scalar (nonlinear) equation for $E^{**}$.

Now $E^{**}$ is known, we can recover all remaining equilibrium components explicitly. Suppose $E^{**}>0$ is a (positive) solution of equation (18) that satisfies $E^{**}<\Lambda /(\sigma +\mu )$. Then compute

$${\lambda }^{**}={\displaystyle \frac{(\sigma +\mu )\mu E^{**}}{\Lambda -(\sigma +\mu )E^{**}}}\mathrm{from}\mathrm{equation}\left(20\right),S^{**}={\displaystyle \frac{\Lambda }{\mu +{\lambda }^{**}}},N^{**}=S^{**}+c{E}^{**}.$$

Then the remaining components are

$$\begin{array}{cc}\hfill I_s^{**}& ={\alpha }_s{E}^{**},\hfill \\ \hfill I_a^{**}& ={\alpha }_a{E}^{**},\hfill \\ \hfill R^{**}& =\rho E^{**}={\displaystyle \frac{\gamma {\alpha }_s+\delta {\alpha }_a}{\mu }}{E}^{**},\hfill \\ \hfill B^{**}& =\zeta E^{**}={\displaystyle \frac{{\omega }_s{\alpha }_s+{\omega }_a{\alpha }_a}{\epsilon }}{E}^{**}.\hfill \end{array}$$

All of these are nonnegative provided $E^{**}\ge 0$. The admissibility condition for this analysis is $0<E^{**}<{\displaystyle \frac{\Lambda }{\sigma +\mu }}$, which ensures the denominator in equation (18) is positive and $S^{**}>0$.

4. Numerical Simulations

Here, we illustrate the theoretical results from the preceding section through numerical simulations of models (1) and (2). Unless specified otherwise, all simulations adopt the parameter values in Table 2.

For numerical simulation of the stochastic AI model, we employed the Euler–Maruyama scheme, which is the standard first–order discretization method for Itô stochastic differential equations [34]. At each time step, the deterministic drift terms of the model were updated using $\Delta t$, while the stochastic perturbations were incorporated through multiplicative noise components of the form ${\sigma }_X{X}_t\Delta W_t$. Let $\Delta t>0$ be the step size and define discrete times

$$t_k=k\Delta t,k=0,1,\dots ,N,$$

where $N=T/\Delta t$ is the total number of steps up to a final simulation time T (500 days). For a Wiener increment we use the discrete approximation

$$\Delta W_k=W\left(t_{k+1}\right)-W\left(t_k\right)\sim \mathcal{N}(0,\Delta t),$$

and in computations we have

$$\Delta W_k=\sqrt{\Delta t}{Z}_k,Z_k\sim \mathcal{N}(0,1).$$

Let $X_k$ denote the numerical approximation of $X\left(t_k\right)$. Euler–Maruyama update Rule applied to equation (4) is

$$X_{k+1}=X_k+f_X\left(X_k\right)\Delta t+{\sigma }_X{X}_k\Delta W_k.$$

The numerical simulation for the DFE analysis shown in Figure 2, compares the deterministic solution of the model (2) with a stochastic realization that incorporates multiplicative environmental and demographic noise. Under the chosen parameter set (Table 2) corresponding to ${\mathcal{R}}_0<1$, the system approaches to disease-free equilibrium, and all infected compartments vanish over time. The deterministic trajectories display smooth exponential decay in the exposed and infectious classes and a monotone return of the susceptible population to its disease-free level. The stochastic realization follows the same qualitative pattern but exhibits sample-path fluctuations. While early dynamics show small irregular oscillations around the deterministic decline in $E,I_a,I_s$ and B, these perturbations do not lead to secondary outbreaks.

Across all infected compartments $(E,I_s,I_a)$ and the environmental viral load $\left(B\right)$, both deterministic and stochastic paths rapidly decline toward zero, confirming the inability of the pathogen to persist within the host population and also in the environment. The small stochastic fluctuations do not enable infection rebound. For susceptible population, despite large sample-path variability, the stochastic path does not show a sustained epidemic wave. It fluctuates around the DFE and remains bounded. Both infectious classes, symptomatic and asymptomatic, decay rapidly in the stochastic simulation, consistent with a negative noise-adjusted growth rate and confirming that random perturbations are insufficient to sustain transmission. Stochastic noise causes small transient departures from the deterministic decline, but the exponential decay dominates. A transient peak exists in the environment early on as infectious hosts shed pathogen into the environment, then decays toward zero. Environmental loading transiently rises but is insufficient to maintain transmission long term due to environmental decay and also low infection prevent sustained environmental forcing. Thus, when $({\mathcal{R}}_0<1)$, neither host-to-host transmission nor environmental transmission can sustain endemic infection in cattle under the current parameter settings.

Figure 3 illustrates how the model behaves under an endemic equilibrium over a 500-day period. All model parameters remain the same as in the baseline disease-free analysis, except for an increase in the asymptomatic transmission rate ${\beta }_a$, which elevates the basic reproduction number to ${\mathcal{R}}_0=3.1945$, pushing the system from a disease-free regime into a stable endemic state. This modification strengthens the indirect and silent spread of infection within the host population, resulting in a higher basic reproduction number. Now the disease no longer dies out; instead, it persists in the population, and the system settles into an endemic state.

In this model, asymptomatic infectious individuals play a crucial role because they continue to interact with susceptible hosts without displaying clinical symptoms. As a result, they remain undetected, untreated, and unrestricted in movement, creating a hidden but highly effective transmission pathway. When ${\beta }_a$ is increased, these individuals contribute more strongly to new infections, amplifying the force of infection. The results shown in Figure 3 also support the condition for $0<E^{**}<{\displaystyle \frac{\Lambda }{\sigma +\mu }}$ in the analysis of endemic equilibrium based on the parameter values as defined in Table 2.

4.1. Sensitivity Analysis

The extinction or persistence of the epidemic depends on the threshold quantity ${\mathcal{R}}_0$. We performed a global sensitivity analysis of ${\mathcal{R}}_0$, as defined in (6), using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) techniques [37,38]. The results are presented in Figure 4, where the bars marked with a red dot indicate the parameters that exert a significant influence on ${\mathcal{R}}_0$. A positive (negative) PRCC value indicates a direct (inverse) relationship with ${\mathcal{R}}_0$, meaning that an increase in the parameter results in an increase (decrease) in ${\mathcal{R}}_0$. For instance, ${\beta }_a,{\beta }_s$, and $\sigma$ have positive values, so increasing any of these parameters increases the value of ${\mathcal{R}}_0$. This suggests that larger values of these parameters contribute to a higher infection burden. Conversely, the parameters $\delta ,\gamma$, and $\nu$ have negative PRCC values, thus, increasing their values results in reduction of ${\mathcal{R}}_0$. Some results of the sensitivity analysis are in simulations in Figure 5 and Figure 6 for two values of ${\beta }_a$ and $\gamma$, respectively.

Figure 5 presents model simulations for two different values of asymptomatic transmission rate ${\beta }_a$. When ${\beta }_a=0.7$, the deterministic trajectory reaches a peak of about 280 symptomatic cattle within roughly 70 days, and the stochastic realizations follow a similarly rapid rise. In contrast, for ${\beta }_a=0.09$, the deterministic peak is lower, at around 65 symptomatic cattle, and occurs later, between 200 and 250 days, with the stochastic trajectories exhibiting greater variability. Consequently, higher transmission leads to earlier and larger outbreaks, whereas lower transmission results in smaller, delayed peaks accompanied by increased stochastic variability.

As illustrated in Figure 6, changing the recovery rate $\gamma$ markedly alters the course of the epidemic. For $\gamma =0.4$, the deterministic trajectory peaks at about day 72 with roughly 75 symptomatic cattle, and the stochastic sample paths remain closely grouped around this relatively small outbreak, indicating rapid infection clearance. In contrast, for $\gamma =0.05$, the deterministic curve increases more gradually but attains a much higher peak of approximately 350 symptomatic cattle around day 85, with the corresponding sample paths displaying greater dispersion. Consequently, a higher recovery rate reduces epidemic magnitude and shortens its duration, whereas a lower recovery rate extends the infectious period, intensifies transmission, and yields larger, later peaks.

4.2. Varying Noise Intensity

Here we examine how changes in the noise intensity, ${\sigma }_k$ (for $k=S,E,I_a,I_s$), affect the epidemic dynamics. The outcomes for three distinct noise intensities are presented in Figure 7. When noise is low (${\sigma }_k=0.010$), the system behaves almost deterministically: the epidemic trajectories for the asymptomatic, symptomatic, and susceptible compartments are smooth, with clearly defined and predictable peaks. The final epidemic size remains tightly clustered around the deterministic value, and the timing of the main outbreak is highly reproducible across realizations. As the noise level increases to ${\sigma }_k=0.050$, moderate stochastic fluctuations introduce visible variability among sample paths. The peak incidence in both asymptomatic and symptomatic groups may be either enhanced or attenuated compared with the deterministic scenario, and the timing of the epidemic peak becomes more variable. Consequently, the distribution of outbreak sizes and durations broadens, complicating both epidemic forecasting and the timely deployment of control strategies.

At high noise levels (${\sigma }_k=0.100$), stochasticity strongly shapes the epidemic course. Sample paths can exhibit multiple irregular peaks, extended low-level circulation, or even premature fade-outs. The symptomatic compartment, in particular, may show large, erratic surges that are not captured by the deterministic model. The susceptible population also displays a much more uncertain trajectory, with some realizations experiencing rapid depletion and others retaining a substantial susceptible pool for long periods. Overall, increasing noise intensity drives the system away from deterministic predictability and toward a regime where epidemic size, peak magnitude, and timing become highly variable and uncertain.

5. Conclusion

This study developed and analysed a stochastic model to explore the transmission dynamics of avian influenza (AI) in cattle populations with both direct and environmental routes of infection. Although AI infections in dairy cattle have only recently gained international attention, especially after the 2024 H5N1 detections in North American herds [14], this model captures the dual role of direct contact transmission and indirect exposure through contaminated environment.

A key insight from the model is the contribution of asymptomatic infections. As seen in livestock influenza systems [21,39], a substantial proportion of infections may be mild or subclinical, yet still infectious. This justifies the separation of symptomatic and asymptomatic infectious states. The model shows that the asymptomatic class can sustain transmission even when symptomatic prevalence is low, particularly when coupled with environmental shedding as seen in Figure 7 (b) and (c). This aligns with recent reports suggesting subclinical AI infections may occur in cattle herds [21], making detection and control more challenging.

As seen in Figure 3, an increase in the value of the asymptomatic transmission rate ${\beta }_a$ leads to the model showing sustained infection and a stable endemic state. The results agree with the sensitivity analysis results in Figure 4. This is biologically feasible for cattle herds, especially in situations with repeated contamination of shared water sources, milking systems, holding yards, or slurry pits.

The introduction of stochasticity into the model plays a crucial role in capturing the inherent variability present in real livestock systems. In contrast to purely deterministic dynamics, which predict smooth trajectories and fixed thresholds, the stochastic model reveals how random fluctuations can substantially alter both the short-term patterns and long-term outcomes of avian influenza transmission in cattle population (see Figure 5). The simulation results show that inclusion of randomness in cattle compartments make the system fluctuates around the endemic equilibrium as seen in Figure 6. This is important because cattle farms may differ in herd size, movement, milking schedule, and hygiene, all of which introduce randomness that deterministic models cannot capture.

Varying the noise intensities across the compartments generally amplifies variability in disease prevalence, but the effects are not uniform across compartments. Instead, the combined action of noise modifies how infection enters, circulates, and persists in the herd. When the overall noise level is low, the stochastic trajectories remain close to the deterministic solution: outbreaks follow predictable curves, and the infection either dies out or settles around a stable endemic level depending on the value of ${\mathcal{R}}_0$. However, as noise intensities increase, the system becomes progressively more irregular, and deviations from the deterministic baseline are amplified, as seen in Figure 7.

One of the most important consequences of higher stochasticity is the increased likelihood of early extinction when the system is near the epidemic threshold. When ${\mathcal{R}}_0>1$ and with varying noise intensities, the infection does not approach a single equilibrium. Instead, the system fluctuates around the endemic level, sometimes producing large transient spikes driven by random increases in viral shedding or environmental contamination. This behaviour is biologically feasible because cattle operations are subject to changes in weather affecting virus survival, variability in shedding between cattle, irregular contact patterns, and differences in cleaning routines; all of these factors contribute to environmental and demographic noise [40].

These results show the importance of stochastic effects when assessing outbreak risks in cattle settings. In real situations, herds rarely exhibit smooth, deterministic dynamics; instead, their disease patterns reflect the cumulative effect of multiple sources of randomness. Incorporating noise into the model results in a more realistic representation of AI dynamics and highlights the importance of measures that reduce environmental variability, improve hygiene, and limit contact patterns in case of an outbreak.

Our findings of with-in herd analysis suggested that transmission rates of asymptomatic and symptomatic cattle, and noise intensity played a significant role in the spread of HPAI disease. This research work is not exhaustive, but it is expected that this would be helpful for livestock farmers, health professionals, and policymakers to take necessary measures against the transmission dynamics of HPAI.

6. Limitations and Future Scope

Although the stochastic model developed in this study provides useful insights into the transmission dynamics of HPAI in cattle population, it has some limitations. The model simplifies heterogeneity between animals by assuming homogeneous mixing, and uniform susceptibility, exposure, and shedding rates. In practice, cattle may differ in age, behaviour, housing conditions, immune responses, and viral shedding patterns. Such heterogeneity can significantly affect transmission, particularly for asymptomatic carriers who may shed virus intermittently or at lower levels. Incorporating age-structured or behaviour-based subpopulations could therefore lead to more accurate predictions. Also, several parameter values were selected from literature on avian influenza in other species due to limited data on cattle-specific AI dynamics. While this approach provides biologically reasonable ranges, it inevitably introduces uncertainty. Direct parameterisation using cattle outbreak data would improve the reliability and predictive accuracy of the model, but such datasets are currently limited.

In future, we will extend this model by incorporating spatially structured systems to capture heterogeneity between farms. Also, we will consider integrating available HPAI longitudinal data to obtain more accurate and context-dependent parameter estimates. Another important extension is evaluating control strategies, such as vaccination and movement restrictions. Stochastic optimal control and cost-effectiveness frameworks could help in understanding avian influenza transmission in cattle and support more resilient disease management strategies.