Probability of a zoonotic spillover with seasonal variation.

Zoonotic infectious diseases are spread from animals to humans. It is estimated that over 60% of human infectious diseases are zoonotic and 75% of them are emerging zoonoses. The majority of emerging zoonotic infectious diseases are caused by viruses including avian influenza, rabies, Ebola, coronaviruses and orthohantaviruses. Spillover of infection from animals to humans depends on a complex transmission pathway, which is influenced by epidemiological and environmental processes. In this investigation, the focus is on direct transmission between animals and humans and the effects of seasonal variations on the transmission and recovery rates. Fluctuations in transmission and recovery, besides being influenced by physiological processes and behaviors of pathogen and host, are driven by seasonal variations in temperature, humidity or rainfall. A new time-nonhomogeneous stochastic process is formulated for infectious disease spread from animals to humans when transmission and recovery rates are time-periodic. A branching process approximation is applied near the disease-free state to predict the probability of the first spillover event from animals to humans. This probability is a periodic function of the time when infection is introduced into the animal population. It is shown that the highest risk of a spillover depends on a combination of animal to human transmission, animal to animal transmission and animal recovery. The results are applied to a stochastic model for avian influenza with spillover from domestic poultry to humans.

Introduction

Zoonotic infectious diseases, transmitted from wild or domestic animal hosts to humans, result in over 60% of human infectious diseases with 75% of them emerging diseases (Guo et al., 2020; Jones et al., 2008; Karesh et al., 2012; Taylor et al., 2001; WHO, 2020b). Many emerging zoonotic diseases are caused by viruses such as avian influenza virus, Ebola virus, rabies virus, severe acute respiratory syndrome coronavirus-1 (SARS-CoV-1), SARS-CoV-2, Middle East respiratory syndrome coronavirus (MERS-CoV), West Nile virus, Nipah virus, Hendra virus and hantaviruses (Han et al., 2015; Poen et al., 2019). A zoonotic spillover occurs when the pathogen from an infected animal host enters a human host, either directly from a natural reservoir, an intermediate animal host, or indirectly from virus in the environment. For example, avian influenza often spills over to humans from an intermediate host (domestic poultry) instead of directly from the natural reservoir (waterfowl including ducks and geese) (Poen et al., 2019; Vandegrift et al., 2010). For hantaviruses, spillover into humans generally occurs via contact with infected excreta from the natural reservoir, such as rats, mice or voles (Jonsson et al., 2010).

Spillovers of public health concern are those involving new pathogens that enter the human population. Animal reservoirs vary considerably, but bats and rodents are the majority of natural reservoirs for viral zoonoses originating in the wild (Calisher et al., 2006; Han et al., 2015; Luis et al., 2013). The natural reservoirs for rabies virus, coronaviruses, Ebola virus, Nipah virus and Hendra virus are bats (Calisher et al., 2006; Cui, Eden, Holmes, & Wang, 2013; Guo et al., 2020; Han et al., 2015; Lau et al., 2005). The Centers for Disease Control and Prevention (CDC), U.S. Department of the Interior and U.S. Department of Agriculture prepared a list of eight zoonotic diseases of most concern to the United States. These eight zoonoses include zoonotic influenza, salmonellosis, West Nile fever, plague, emerging coronavirus diseases, rabies, brucellosis and Lyme disease (CDC, 2019a). Four of the eight pathogens are viruses.

Zoonotic diseases can be transmitted to humans directly or indirectly. Direct or indirect transmission requires contact between the animal source and humans or contact with the pathogen in the environment. These contacts are often driven by environmental factors, especially seasonal variations in temperature, humidity and rainfall that impact food resources and habitat e.g., (Glass et al., 2000; Karesh et al., 2012; Martinez, 2018; Plowright et al., 2017).

With the exception of vector-host epidemic models, much of the modeling efforts for emerging viral zoonoses have concentrated on the disease dynamics within a single population, either humans or the natural reservoir. There are a number of zoonotic modeling studies that include at least two species, such as models for avian influenza (Gumel, 2009; Iwami et al., 2007; Liu et al., 2017; Royce and Fu, 2020; Tuncer & Martcheva, 2013; Vaidya & Wahl, 2015; Zhang et al., 2019), rabies (Asamoah, Oduro, Bonyah, & Seidu, 2017; Huang et al., 2019; Ruan, 2017), hantaviruses (Allen et al., 2009; Sauvage et al., 2007) and COVID-19 (Chen et al., 2020). Most of these models are systems of ordinary differential equations (ODEs). A few stochastic zoonotic disease models include human spillover or seasonal variability, e.g., (Breban, Drake, Stallknecht, & Rohani, 2009; Keeling & Gilligan, 2000; Singh et al., 2014; Vaidya & Wahl, 2015; Voinson et al., 2018). Breban et al. (2009) and Vaidya and Wahl (2015) included seasonal variability in models for avian influenza in migratory waterfowl, as they move between breeding and wintering sites. Keeling and Gilligan (2000) modeled spillover of bubonic plague for a metapopulation of rats, fleas and humans that included seasonal forcing in the flea carrying capacity. Stochastic simulations of the metapopulation model showed sporadic human outbreaks from subpopulations of endemic rat reservoirs. General stochastic models for spillover to humans from an animal reservoir were studied by Voinson et al. (2018) and Singh et al. (2014). These investigations applied Markov chain and branching process theory to susceptible-infected-recovered (SIR) models to study the number and size of outbreaks. Singh et al. (2014) also investigated the dynamics during the early epidemic stages through estimation of the probability of and time to first human spillover. We apply similar techniques, as in Singh et al. (2014) and Voinson et al. (2018), to investigate the early epidemic dynamics. In addition, we extend their models by including seasonal variability in the parameters and by considering more complex settings.

In this investigation, we develop mathematical methods for study of the first spillover event to humans via direct transmission either from a natural reservoir or an intermediate host. We focus on viral zoonoses. However, the methods can be extended to other zoonotic pathogens and to indirect exposure. The following section describes the underlying SIR ODEs from which the transition rates for the time-nonhomogeneous Markov chain are defined. The branching process approximation is described in Section 3. Then, in Section 4, these new mathematical methods are applied. Numerical examples illustrate how the periodic parameters affect the probability of and time to spillover as a function of the time when infection is introduced into the animal population. In Section 5, the methods are applied to several avian influenza models with spillover from domestic poultry to humans (Tuncer & Martcheva, 2013), where the SIR framework is also generalized to include a latent stage. In the last section, the results are summarized and the public health implications are discussed.

Spillover model

A system of SIR ODEs serves as a framework for formulation of a stochastic model for spillover. Since we are considering the first spillover event, the human disease dynamics after the spillover are not modeled. Whether an outbreak occurs in the human population depends on human-to-human transmission, human recovery rate (Singh et al., 2014; Voinson et al., 2018) and seasonal effects (Nipa & Allen, 2020).

Deterministic spillover model

A compartmental diagram of the deterministic spillover model is given in Fig. 1. The animal population includes susceptible, infectious and recovered animals, denoted as S, I and R, respectively. Similarly, S, I and R denote susceptible, infectious and recovered humans, respectively. The parameter β(t) is the transmission rate from animal to animal, γ(t) is the recovery rate of an infected animal and β(t) is the transmission rate of animal to human. Each of these rates are time-periodic with a common period T. The population sizes are constant, N = S + I + R and N = S + I + R. For the human population, we do not consider the recovered state, because our interest is in the first spillover event in the human population. The ODE model for the animal and the human populations is described by the following system:

Fig. 1

Compartmental diagram of the SIR ODE model for a spillover event from an animal host (subscript a) to a human host (subscript h).

We assume the parameters are continuous and periodic for t ∈ (−∞, ∞) with common period T > 0:one or more of the rates may be constant or have a fundamental period other than T.

Stochastic spillover model

The susceptible, infectious and recovered animal and human variables are discrete random variables satisfyingthe values of the random variables are denoted with lower case letters: s, i, r, s and i. The initial conditions are at a fixed time t0 ∈ [0, T]. They are I(t0) > 0 and R(t0) = I(t0) = R(t0) = 0 with S(t0) = N − I(t0) and S(t0) = N. To record the first spillover event, the stochastic process continues until a time t (t > t0) such that either I(t) = 0 or I(t) = 1. We compute the probability that I(t) = 1 before I(t) = 0.

The transition probabilities that lead up to a spillover are described in Table 1. Given I(t0) > 0 and I(t0) = 0, three possible events may occur: (i) an infected animal infects another animal, (ii) an infected animal recovers or (iii) an infected animal infects a human. Event (iii) is a spillover. The sum of the transition rates is

Table 1

Transition probabilities for a zoonotic spillover, r(t)Δt + o(Δt).

Description	Event	(ΔIa, ΔIh)	Probabilities
Animal infection	(Sa, Ia) → (Sa − 1, Ia + 1)	(1, 0)	β_aa(t)i_as_aN_aΔt+o(Δt)
Animal recovery	(Ia, Ra) → (Ia − 1, Ra + 1)	(−1, 0)	γ(t)iaΔt + o(Δt)
Human infection	(Sh, Ih) → (Sh − 1, Ih + 1)	(0, 1)	β_ah(t)i_as_hN_hΔt+o(Δt)
No change	(Sa, Ia, Sh, Ih) → (Sa, Ia, Sh, Ih)	(0,0)	1 −Σ(t)Δt + o(Δt)

To numerically approximate the Markov time-nonhomogeneous process, we use a Monte Carlo approximation based on Table 1. We consider a discrete set of time points with t = nΔt for sufficiently small time step Δt so that the probabilities r(t)Δt, i = 1, 2, 3 and 1 −Σ(t)Δt lie in [0, 1] for t > 0. A uniform random number is used to select the event, if any, that occurs at time t+1. An alternate method based on a Gillespie-type algorithm is discussed in Appendix A.

In Fig. 2, five sample paths of the infected animal population are graphed for two different initial conditions. The two initial conditions are I(t0) = 1 for t0 = 0, 1 ∈ [0, T], where the period length is T = 4. The sample paths are graphed until the time t, when either I(t) = 0 or I(t) = 1, t > t0, whichever event occurs first. As can be seen in Fig. 2, there are significant differences in the occurrence of a spillover. At t0 = 0, transmission rates β(0) and β(0) and recovery time 1/γ(0) are at their peak values; four out of five sample paths result in a spillover. At t0 = 1, the values of β(1), β(1) and 1/γ(1) have much smaller values than at t0 = 0 and only one of the five sample paths results in a spillover.

Fig. 2

Five sample paths are graphed for the infected animal population in the stochastic spillover model until either I(t) = 0 or I(t) = 1, t > t0. Parameter values are N = 1000 = N, , i = aa, ah, , , and and initial conditions S(t0) = 999, I(t0) = 1, S(t0) = 1000 and I(t0) = R(t0) = R(t0) = 0 for t0 = 0 or t0 = 1.

Branching process approximation

We apply a Markov branching process approximation to the infected animal population, variable I, at the disease-free equilibrium (DFE), S = N and S = N. In particular, we assume I follows a birth-death-killing process, as defined by Karlin and Tavaré (1982), with state space {0, 1, ⋯ }∪{K}, where state K is the first human spillover event. The two states, I = 0 and I = K (I = 1), are the only absorbing states for I. The process begins at a fixed time t0 and the process stops at time t > t0, if either I(t) = 0 or I(t) = 1. Karlin and Tavaré (1982) assumed the birth-death-killing process was time-homogeneous with constant parameters. The time-homogeneous process has been applied by Singh et al. (2014) to spillover events for SIR epidemic processes in humans and animals. We extend their work to a time-nonhomogeneous process for periodic transmission and recovery rates and to more general epidemic settings. We summarize some of the theory and describe how the probability of spillover is computed.

The assumptions in the Markov branching process for I are that the population size is unbounded, I ∈{0, 1, 2, …}, the infected animals at time t are independent of those infected prior to time t (Markov property) and infected animals produced after a time t from infected animals at time t are independent of each other (Harris, 1963). To ensure the branching process is a good approximation of the time-nonhomogeneous process at the human-animal interface, the value of N must be sufficiently large and the initial number of infected animals must be sufficiently small so that the process is near the DFE.

For the Markov branching process, the events in Table 1 simplify to those given in Table 2. Given I(t0) = i, the branching process is used to estimate the probability of animal disease extinction before the first human spillover. We denote this probability as . The probability of the first human spillover event is . The derivation of this expression is described below.

Table 2

Transition probabilities for a branching process approximation of spillover.

Event	(ΔIa, ΔIh)	Probabilities
(i)	(1, 0)	βaa(t)iaΔt + o(Δt)
(ii)	(−1, 0)	γ(t)iaΔt + o(Δt)
(iii)	(0, 1)	βah(t)iaΔt + o(Δt)

The expression for the probability of the first human spillover comes from the more general multitype branching process at the human-animal interface (Singh et al., 2014). The generating function for the interface involves both I and I. Given the transition probability,the probability generating function for (I, I) equals

As our focus is the first human spillover, we restrict to the case I(s) = 0 for s ∈ [t0, t], i.e., j = l = 0 in equations (5), (6). With these restrictions, the simplified expressions for the transition probability and the generating function are defined as follows:in particular,

Applying the probabilities in Table 2, we derive the backward Kolmogorov differential equations for the restricted process. It follows that the transition probability for a small change in the initial time Δt0 > 0 equalssubtracting p,(t0, t), dividing by Δt0 and letting Δt0 → 0+ yieldthe probability of animal disease extinction at time t (with no human spillover) from one infected animal at time t0 equals

It follows from the branching process assumption of independence that (Harris, 1963). This property also applies to p,0(t0, t) when u = 0. In particular,using the identity (9) in the right side of equation (7) when i = 1 and k = 0, setting p0,0(t0, t) = 1 and simplifying yieldwhere f is defined asthe generating function f has the term β(t0) in the denominator implying animal disease extinction is conditional on no human spillover.

The partial differential equation in (10) does not depend explicitly on t. Therefore, fix t and make a change of variable s = t − t0 in equation (10) so that P(s) = p1,0(t0, t). The following equation is solved backward in time:application of Theorem B.1 and Proposition B.1 in Appendix B show that P(s) converges monotonically and uniformly to a periodic solution Φ(s) on s ∈ [0, T]. In particular,

For this branching process, the solutions of the backward and forward Kolmogorov differential equations for I(t0) = 1 are in agreement (Feller, 1940). After a change of variable, the asymptotic solution in (13) can be applied forward in time to p1,0(t0, t) as t → ∞:thus, given I(t0) = 1, the asymptotic periodic probability for disease extinction in the animal population equals

It follows that the asymptotic probability of spillover for I(t0) = 1 equals . Finally, the branching process assumption (9) leads to the expression for the asymptotic probability of spillover when I(t0) = i:

To compute numerically, first we solve equation (12) for sufficiently large s > 0 to obtain an accurate approximation of the periodic solution Φ(s). Second, we define and finally, we apply formula (15) which gives the desired probability of spillover.

In some special cases, the probability of spillover may have a trivial periodic solution, i.e., a constant solution. One obvious case is when the coefficients β, β and γ are constant. Another case is when f(u, t0) in (11) is independent of t0, f(u, t0) ≡ f(u). More generally, a trivial periodic probability of spillover occurs if there exists q ∈ (0, 1) such that f(q, t0) = q for all t0 ∈ [0, T]. In these cases, the probability of spillover is . These results are stated in the next theorem. The Proof is given in Appendix B.

Assume the transmission and recovery rates are positive, continuous and periodic functions defined on with common period T > 0. In addition, assume that there exists a constant q ∈ (0, 1) such that f(q, t) = q for all t ∈ [0, T]. Then q is the unique fixed point with this property. In addition, given I(t0) = i, the asymptotic probability of spillover from the branching process approximation is constant,for all t ∈ [0, T]. In the particular case that the following ratios are constant,an explicit expression for q is given by

When the parameter values β, β and γ are constant, the two ratios and can be characterized as reproduction numbers. The ratio is the basic reproduction number for the animal population, i.e., the average number of secondary infections caused by one infected animal during the animal’s infectious period (van den Driessche & Watmough, 2002). The ratio is the spillover reproduction number, i.e., the average number of spillover infections caused by one infected animal during the animal’s infectious period.

Numerical examples

Several numerical examples are presented to illustrate how the probability of spillover depends on the three parameters, β, β and γ. In the first example, we illustrate this dependence when the parameter values are constant and therefore, the two reproduction numbers, and , in Theorem 3.1 are constant. In Fig. 3, the probability of spillover is graphed as a function of these two reproduction numbers. An increase in either of these two reproduction numbers increases the probability of spillover. Interestingly, even if the basic reproduction number for the animal population is less than one, (no major outbreaks), there may still be a spillover in the human population. But in this latter case, the spillover reproduction number must be sufficiently large. When , is sensitive to the value of . When , there are either major or minor outbreaks. For I(t0) = 1, the probability of a minor animal outbreak is while a major outbreak is (Whittle, 1955). Hence, the probability of no spillover, conditional on a minor animal outbreak, is and the probability of a spillover, conditional on a minor animal outbreak, is . Fig. 3 illustrates the difference between the probability of a spillover when there is a minor or a major outbreak. Fig. 3 agrees with Fig. 2 in Singh et al. (2014).

Fig. 3

The probability of spillover, , as a function of and when parameter values are constant (Theorem 3.1).

In the remaining numerical examples, we assume periodic transmission and recovery rates and compare the estimate of the probability of spillover in equation (15) to the value obtained from the sample paths of the time-nonhomogeneous process. For illustration purposes, suppose the transmission and recovery rates are trigonometric functions with period T = 4 seasons (one season ≈ 90 days):and other parameter values are given in Table 3. The rates have units per season. The average values of the transmission and recovery rates are , i = aa, ah and , respectively. The absolute values of the ε, i = aa, ah, g are the amplitudes of the seasonal rates.

Table 3

Parameter values for the stochastic spillover model.

Parameter	Value	Parameter	Value
Na	1000 animals	εi, i = aa, ah	±0.9
Nh	1000 humans	β¯_aa	4/season
T	4 seasons	β¯_ah	1/season
εg	±0.9	γ¯_a	6/season

In the numerical simulations of the Markov time-nonhomogeneous process, the events are defined as in Table 1. A total of 104 sample paths are generated for each of five different initial time points t0 = 0, 1, 2, 3, 4, during the period T = 4. The proportions of the 104 sample paths that result in a human spillover provide numerical estimates of the probabilities of a spillover. The remaining proportions are those sample paths that result in disease extinction in the animal population with no human spillover.

In Fig. 4, it can be seen that the estimate of the spillover probabilities from the branching process solution of equation (15) are in good agreement with the numerical simulations of the time-nonhomogeneous process (blue circles) when I(t0) = 1 or 5. The peak values of the probabilities of spillover are close to, but shifted left of the peak values of one of the transmission rates, either β(t0) (rows (b) and (c)) or β(t0) (row (d)).

Fig. 4

Probabilities of spillover for different combinations of the periodic parameters γ(t), β(t) and β(t). Graphs of the periodic parameters are in the left panels. The branching process estimate of probability of spillover is graphed in the middle and right panels when either one or five infected animals are introduced, I(t0) = 1 or 5 (black curves). Probabilities of spillover from simulation of 104 sample paths of the time-nonhomogeneous process are marked by blue circles at t0 = 0, 1, 2, 3, 4. Parameter values are in Table 3 and periodic parameters are in equations (18), (19). (a) First row: ε = ε = ε = 0.9. (b) Second row: ε = ε = 0.9 and ε = −0.9. (c) Third row: ε = −0.9 and ε = ε = 0.9. (d) Fourth row: ε = ε = −0.9 and ε = 0.9.

In Fig. 5 are the graphs of the mean and standard deviation for the time until the first spillover when I(t0) = 1 or 5, corresponding to each of the four cases in Fig. 4. The times at which infection is initiated are t0 = 0, 0.5, 1, …, 3.5. The longest mean times and greatest variability occur in case (a) when the three periodic functions have the same shape and they are all close to zero (slow change in the dynamics). The shortest mean times occur for case (c) when the transmission and the recovery rates have the reverse amplitudes, i.e., the recovery rate is fast and the transmission rates are slow.

Fig. 5

Plots of the mean and standard deviation for the time until the first spillover event for the transmission and recovery rates given in Fig. 4 (a)–(d) when I(t0) = 1 or 5. The times when infected animals are introduced are t0 = 0, 0.5, 1, …3.5.

The average probability of a spillover is computed as follows:the average probabilities for I(t0) = 1 or 5 are computed for each of the periodic spillover probabilities graphed in Fig. 4. They are summarized in Table 4. The results depend on the values of the periodic parameters and the initial number of animals infected. For the particular examples in Fig. 4, the average probabilities of spillover with periodic rates in (b), (c) or (d) may be greater or less than the values in case (a) which correspond to constant rates.

Table 4

The average probability of spillover when the initial infected animal population is I(t0) = 1 or 5, based on the periodic spillover probabilities graphed in Fig. 4 (a)–(d).

Fig. 4	P¯_spill(1)	P¯_spill(5)
(a)	0.250	0.763
(b)	0.397	0.785
(c)	0.373	0.626
(d)	0.304	0.676

The role of the transmission and recovery rates in predicting the peak values of the spillover probabilities is investigated. We let the recovery rate be constant and assume the transmission rates are periodic. In Fig. 6, the probabilities of spillover are graphed for four different constant recovery rates γ = 2, 4, 6, 8. The transmission rates, given in equation (18), are periodic with fundamental period T = 4. In Fig. 6 (a), the transmission rates are synchronized in time with β(t) = 4β(t) (similar to Fig. 4 (a) and (c)). In Fig. 6 (b), the transmission rates are desynchronized in time, with β(t) shifted right T/2 units of β(t) (similar to Fig. 4 (b) and (d)). The four graphs, top to bottom, in the right panels of Fig. 6 correspond to the four recovery rates: γ = 2, 4, 6, 8. As the parameter γ is increased, the minima and maxima of the probability of spillover switch from being dominated by the extrema of β(t) to those of β(t). When γ increases, the average reproduction numbers decrease: and . In Fig. 6, if , the times of the peak values of the animal to animal transmission rate are more closely related to the times of highest risk of spillover than the time of the peak values of the animal to human transmission, whereas if , the reverse is true. For and , the spillover probability is constant, (Theorem 3.1 and Figure B.1 in Appendix B). A similar phenomenon can be observed if γ = 6 is fixed and the mean value of the animal to human transmission rate, , is increased from values less than one to greater than one.

Fig. 6

Probability of spillover when γ is a fixed constant and β(t) and β(t) are periodic. The periodic transmission rates are graphed in the left panels and the corresponding probabilities of spillover, as computed from the branching process approximation, are graphed in right panels for four different values of γ = 2, 4, 6, 8. In the two panels on the right, the times at which there is greatest risk of a spillover (maximum values of ) are marked with blue plus sign + and the times of lowest risk (minimum values of ) are marked with a red dot ⋅.

Application to avian influenza

Avian Influenza (AI) is caused by influenza A viruses that are classified according to two surface proteins, hemagglutinin (H, 16 subtypes) and neuraminidase (N, 9 subtypes) (CDC, 2017). The AI viruses that infect humans are generally associated with either H5, H7 or H9 subtypes (CDC, 2019b). Outbreaks around the world of H5N1 from 2003 to July 2020, recorded by the World Health Organization, resulted in 861 confirmed human cases and 455 deaths (WHO, 2020a). An H7N9 outbreak in China from 2013 to 2015, resulted in 568 confirmed human cases and 212 deaths (WHO, 2015). A few sporadic and mild human cases have been reported for H9N2 (CDC, 2019b). Currently, none of these AI viruses have sustained human-to-human transmission (WHO, 2015; 2020a). The natural reservoir for AI is wild birds, specifically migrating waterfowl. Spread to humans often comes from an intermediate host, such as domestic poultry. AI viruses are also classified according to low or highly pathogenic AI (LPAI or HPAI) which refers to the severity of disease in domestic poultry (CDC, 2017). Human infection with H5N1 is often associated with HPAI while the other two, H7N9 and H9N2, are LPAI (CDC, 2017). Seasonal outbreaks of HPAI are frequently seen in domestic poultry (OIE, 2018; Tuncer & Martcheva, 2013), especially for HPAI, where seasonal migration patterns of wild birds have a strong connection with outbreaks in domestic poultry (Bui et al., 2016).

A variety of models for AI include transmission between birds and humans but only a few include seasonal variability, e.g., (Liu et al., 2017; Tuncer & Martcheva, 2013; Zhang et al., 2019). Tuncer and Martcheva (2013) applied an ODE model for spread of HPAI H5N1 between domestic poultry and humans and fit seven different models to cumulative number of human cases. The models included seasonality in domestic bird-to-bird transmission, wild bird migration, environmental transmission or a combination of all three types of seasonality. We apply the seasonal HPAI H5N1 model formulated by Tuncer and Martcheva (2013) that provided the best fit to cumulative number of human cases, the model with seasonal bird-to-bird transmission. The model also assumes replacement, removal and disease-related mortality of domestic birds but no recovery. Parameter Λ is the bird replacement rate, μ is the removal rate and ν is the disease-related mortality. In the absence of disease, the bird population is maintained at a constant population size of S = Λ/μ = N. The total human population is N. The seasonal bird-to-bird transmission is β(t). We also consider a case with seasonal bird-to-human transmission, β(t). The ODE model (1) is replaced by the following system:and the human population is modeled by system (2).

A Markov time-nonhomogeneous stochastic process can be defined based on the rates in the ODE models (2)–(21) as in Table 1. In this new stochastic model, there are events for replacement, removal and disease-related mortality in the domestic bird population. Approximation of the stochastic process near the DFE, S = N and S = N, leads to a branching process approximation of the time-nonhomogeneous process. Similar to the derivation of the differential equation (12), assuming I(t0) = 1 and defining s = t − t0, leads to a new differential equation for the probability of no spillover for the branching process approximation of the stochastic model:where the new generating function f equals:

The periodic probability of spillover for I(t0) = 1 is found by solving the differential equation (22) with the generating function in equation (23). To test the branching process estimate for the probability of spillover, we use the parameter values in Tuncer and Martcheva (2013), but we adjust the bird and human population sizes N and N so that they correspond to population sizes on a single large poultry farm. We also adjust so that the average reproduction numbers correspond to those in Tuncer and Martcheva (2013) (Table 5). In our model, population sizes for birds and humans are N = 20, 000 and N = 200 and the average reproduction numbers equal

Table 5

Parameter values for domestic birds and humans. Parameter values are taken from Tuncer and Martcheva (2013) with the exception of N, N and .

Parameter	Poultry	Parameter	Human
Na	20,000 animals	Nh	200 humans
β¯_aa	0.1075/day	β¯_ah	0.1235/day
ν	0.1/day	T	365 days
μ	0.5/T/day

Three different sets of periodic transmission rates for HPAI H5N1 are applied to the domestic bird-human stochastic model:case (i) is applied in Tuncer and Martcheva (2013). The periodic probability of spillover for these three cases are graphed in Fig. 7. There is good agreement between the branching process approximation for t0 ∈ [0, 360] and the numerical simulations of the time-nonhomogeneous stochastic process for initial times t0 = 0, 30, 60, 90, …, 360 when I(t0) = 1 or 2 (104 sample paths at each initial time point). The differences in amplitude and shape of the periodic transmission rates over time have a large impact on the probability of spillover. The average values of the probability of spillover are summarized in Table 6. In case (i), the average values are close to those in the model with constant transmission rates, e.g., and . In cases (ii) and (iii) with more variability in the seasonality, the average values are smaller than in case (i).

Fig. 7

Probability of spillover from the branching process approximation of the time-nonhomogeneous process derived from the ODE model (2)–(21). The transmission rates β and β are graphed in the left panels for each of the three cases (i), (ii) and (iii), equations (24), (25), (26), respectively. Other parameter values are given in Table 5. In the two panels on the right are the probability of human spillover, computed from the branching process approximation when either 1 or 2 infected birds are introduced into the population (black curves). From simulation of the time-nonhomogeneous process at t0 = 0, 30, 60, 90, …, 360, the proportions of sample paths out of a total of 104 that result in a human spillover before death or removal of birds are also plotted (blue circles).

Table 6

Average probability of spillover when the initial domestic bird population equals I(t0) = i, i = 1, 2 based on the branching process approximation with seasonal transmission rates as in cases (i), (ii), (iii).

Case	P¯_spill(1)	P¯_spill(2)
(i)	0.657	0.882
(ii)	0.651	0.875
(iii)	0.564	0.758

Generalization to multiple stages

The results apply to more general spillover models that include multiple stages and seasonality. We illustrate how the spillover probabilities change with inclusion of a latent stage in the animal population. If deaths occur during the latent stage, then the probabilities of spillover will be decreased. Let the latent stage be denoted as E. Let the mortality rate (natural plus disease-related) be μ + α and the transition rate to the infectious stage be δ. The ODE model for the domestic bird population is

The rates in the domestic bird ODE model (27), coupled with the rates in the human ODE model (2), can be used to define transition probabilities for the time-nonhomogeneous stochastic model. The branching process approximation of the stochastic model at the DFE, S = Λ/μ = N and S = N, and the backward Kolmogorov differential equations can be used to derive two differential equations, one for each of the variables E and I (similar to the method for derivation of (10)). A change of variable in the two differential equations, s = t − t0, leads to the following system of differential equations:with initial conditions P1(0) = P2(0) = 0 and the two generating functions

The generating function f1, given E(t0) = 1, accounts for two events. Either E transitions to stage I with probability δ/(δ + μ + α) or dies with probability (μ + α)/(δ + μ + α). The generating function f2, given I(t0) = 1, also accounts for two events, conditional on no spillover. Either I infects another animal resulting in a new exposed animal E with probability β(t0)/(β(t0) + μ + ν + β(t0)) or dies with probability (μ + ν)/(β(t0) + μ + ν + β(t0)). Numerical solution of system (28) leads to an estimate for the periodic probability of spillover:for i = 1, 2, as k → ∞. Thus, , and

In Figure (8) are graphs of the branching process estimates of the spillover probabilities when the transmission rates are given by case (ii). Spillover is initiated by domestic birds with either one exposed bird or one infected bird: or (0, 1). The new parameters are δ = 0.2/day and α = 0.1/day with other parameter values in Table 5.

Fig. 8

Probability of spillover for the AI time-nonhomogeneous stochastic process corresponding to (2)–(27) with transmission rates β and β for case (ii). Parameters are δ = 0.2/day, α = 0.1/day with other parameter values in Table 5. Numerical simulations from the time-nonhomogeneous process show the proportions of sample paths out of 104 that result in a human spillover before death or removal of all of the domestic birds (blue circles) at t0 = 0, 30, 60, 90, …, 360.

The average value of the periodic probability of spillover depends on the initial conditions (E(t0), I(t0)) = (j, i), t0 ∈ [0, T]. It can be computed similar to formula (22) asTable 7 is a summary of the average values for the probability of spillover for cases (i), (ii) and (iii) with either 1 or 2 exposed or infected birds.

Table 7

Average probability of spillover from the branching process approximation of the time-nonhomogeneous process corresponding to ODE model (2)–(27) for cases (i)-(iii) for either 1 or 2 exposed or infected birds introduced into the domestic bird population.

Case	P¯_spill(1,0)	P¯_spill(2,0)	P¯_spill(0,1)	P¯_spill(0,2)
(i)	0.414	0.656	0.624	0.858
(ii)	0.412	0.654	0.621	0.855
(iii)	0.358	0.566	0.539	0.739

The average values of the probabilities of spillover with an exposed stage differ from those without an exposed stage, Table 6 versus Table 7. As the exposed stage includes deaths, average values are smaller if spillover is initiated by exposed rather than by infected birds. Other generalizations of the branching process method are possible, such as spillover from indirect transmission through the environment or transmission involving several species, a natural reservoir, an intermediate host and humans.

Discussion

Seasonality plays a crucial role in the timing of first human spillover in zoonotic diseases (Glass et al., 2000; Plowright et al., 2017; Schmidt et al., 2017). In this investigation, we generalized the mathematical results of Singh et al. (2014) by including the effects of seasonality during the early stages of an epidemic. We applied branching process theory and the backward Kolmogorov differential equations to derive an analytical approximation for the probability of the first human spillover when parameters such as the transmission and recovery rates vary seasonally. The probability of spillover depends on the time during the season when infection is introduced into the animal population. Prediction of the time of highest risk of human spillover depends on how seasonality affects the parameter rates and generally, may not be obtained by observation of a single parameter (Fig. 4, Fig. 6).

The assumption of unbounded population sizes in the branching process limits the application of these results to large population sizes. The branching process estimates should be tested against the Markov time-nonhomogeneous process.

From a public health perspective, prevention and management of spillover require medical and veterinary approaches such as vaccination and treatment as well as ecological considerations that target transmission pathways at the human-animal interface (Karesh et al., 2012; Lloyd-Smith et al., 2009; Sokolow et al., 2019). Seasonality and changing climate are important drivers of this spillover. Understanding the role of seasonality in spillover prevention and management requires models coupled with data. Important questions for models to address are the effects of medical, veterinary and ecological interventions on public health outcomes when transmission is seasonal (Karesh et al., 2012; Lloyd-Smith et al., 2009; Schmidt et al., 2017; Sokolow et al., 2019).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.