key: cord-337992-g4bsul8u
authors: Voinson, Marina; Alvergne, Alexandra; Billiard, Sylvain; Smadi, Charline
title: Stochastic dynamics of an epidemic with recurrent spillovers from an endemic reservoir
date: 2018-11-14
journal: Journal of Theoretical Biology
DOI: 10.1016/j.jtbi.2018.08.017
sha: 
doc_id: 337992
cord_uid: g4bsul8u

Abstract Most emerging human infectious diseases have an animal origin. While zoonotic diseases originate from a reservoir, most theoretical studies have principally focused on single-host processes, either exclusively humans or exclusively animals, without considering the importance of animal to human transmission (i.e. spillover transmission) for understanding the dynamics of emerging infectious diseases. Here we aim to investigate the importance of spillover transmission for explaining the number and the size of outbreaks. We propose a simple continuous time stochastic Susceptible-Infected-Recovered model with a recurrent infection of an incidental host from a reservoir (e.g. humans by a zoonotic species), considering two modes of transmission, (1) animal-to-human and (2) human-to-human. The model assumes that (i) epidemiological processes are faster than other processes such as demographics or pathogen evolution and that (ii) an epidemic occurs until there are no susceptible individuals left. The results show that during an epidemic, even when the pathogens are barely contagious, multiple outbreaks are observed due to spillover transmission. Overall, the findings demonstrate that the only consideration of direct transmission between individuals is not sufficient to explain the dynamics of zoonotic pathogens in an incidental host.

Recent decades have seen a surge of emerging infectious diseases (EIDs), with up to forty new diseases recorded since the 1970s ( Jones et al., 2008 ) . Sixty percent of emerging human infectious diseases are zoonotic, i.e. are caused by pathogens that have an animal origin ( Jones et al., 2008; Taylor et al., 2001 ) . The World Health Organization defines zoonotic pathogens as "pathogens that are naturally transmitted to humans via vertebrate animals". The epidemics caused by EIDs impact the societal and economical equilibria of countries by causing unexpected deaths and illnesses thereby increasing the need for health care infrastructures and by interfering with travel ( Morens and Fauci, 2013 ) . Moreover, the risk of EIDs being transmitted to humans from wildlife is increasing because of the recent growth and geographic expansion of human populations, climate change and deforestation, which all in-crease the number of contacts between humans and potential new pathogens ( Jones et al., 2008; Keesing et al., 2010; Murray and Daszak, 2013 ) . Given most EIDs have an animal origin, it is crucially important to understand how infections spread from animal to human populations, i.e. by spillover transmission.

There is numerous empirical evidence that the epidemiological dynamics of infectious diseases is highly dependent on the transmission from the reservoir (the reservoir will be defined following Ashford's definition (1997) , i.e. a pathogen is persistent in the environment of the incidental host, see Table 1 for details). The start of an outbreak is promoted by a primary contact between the reservoir and the incidental host (i.e. host that becomes infected but is not part of the reservoir) leading to the potential transmission of the infection to the host population. Moreover, multiple outbreaks are commonly observed during an epidemic of zoonotic pathogens in human populations, for instance in the case of the epidemic of the Nipah Virus between 2001 and 2007 ( Luby et al., 2009 ) . With regards to the Ebola virus, some twenty outbreaks have been recorded since the discovery of the virus in 1976 ( De La Vega et al., 2015 ) . This number of outbreaks undoubtedly underestimates the https://doi.org/10.1016/j.jtbi.2018.08.017 0022-5193/© 2018 Elsevier Ltd. All rights reserved.

Definitions of a reservoir from the literature. The reservoir is mostly used as defined by the Centre for Disease Control and prevention (CDC) . Two other definitions have been proposed to clarify and complete the notion of reservoir in the case of zoonotic pathogens. On the one hand, Haydon et al. (2002) define the reservoir from a practical point of view in order to take into account all hosts epidemiologically connected to the host of interest (i.e. target host), to implement better control strategies. On the other hand, Ashford (1997) establishes a more generalizable definition: for a given pathogen there is a single reservoir.

Authors Refs "any animal, person, plant, soil, substance or combination of any of these in which the infectious agent normally lives"

Centre for Disease Control and Prevention (2018) "all hosts, incidental or not, that contribute to the transmission to the target host (i.e. the population of interest) in which the pathogen can be permanently maintained" Haydon et al. (2002) "an ecologic system in which an infectious agent survives indefinitely" Ashford (1997 Ashford ( , 2003 total number of emergences because not all emergences necessarily lead to the spread of the infection from an animal reservoir to the host population ( Jezek et al., 1999 ) . While the reservoir has an important role for causing the emergence of outbreaks, the role of spillover transmission on the incidental epidemiological dynamics is rarely discussed. Empirically, it is generally difficult to distinguish between direct transmission and transmission from the reservoir. Only in the case of non-communicable diseases it is easily possible to measure the importance of the recurrent transmission from the reservoir, since all infected individuals originate from a contact with the reservoir. For instance, the H7N9 virus, for which most human cases are due to a contact with an infected poultry, approximately 132 spillovers have been listed during the epidemic of 2013 ( Zhou et al., 2013 ) . For pathogens that are able to propagate from one individual to another, the origin of the infection can be established according to patterns of contacts during the incubation period ( Chowell et al., 2014; Luby et al., 2009 ). Most often, if an infected individual has been in contact with another infected individual in his recent past, direct transmission is considered as the likeliest origin of the infection. However, both individuals might have shared the same environment and thus might have been independently infected by the reservoir. This leads to overestimating the proportion of cases that result from person-to-person transmission. Moreover, when the pathogen infects an individual and the latter does not produce secondary cases then the detection of emergence is unlikely.

Pathogen spillover is often neglected in epidemiological theoretical models. It is generally assumed that the epidemiological dynamics of outbreaks is driven by the ability of the pathogen to propagate within hosts. For instance, a classification scheme for pathogens has been proposed by Wolfe et al. (2007) , including five evolutionary stages in which the pathogen may evolve ranging from an exclusive animal infection (Stage I) to an exclusive human infection (Stage V) ( Fig. 1 ) ( Wolfe et al., 2007 ) . The intermediate stages are those found for the zoonotic pathogens (Stages II-IV). Lloyd-Smith et al. (2009) , proposed to enhance the classification scheme by differentiating the Stages II-IV by the ability of pathogens to propagate between individuals in the incidental host (i.e. as a function of the basic reproductive ratio R 0 ): the noncontagious pathogens ( R 0 = 0 , Stage II), pathogens barely contagious inducing stuttering chains of transmission (0 < R 0 < 1, Stage III) and contagious pathogens inducing large outbreaks ( R 0 > 1, Stage IV) ( Lloyd-Smith et al., 2009 ) . However, the role of the reservoir is not clearly defined, and spillover effects on the epidemiological dynamics are not discussed.

Only a few models have investigated the dynamics of EIDs by taking into account explicitly the transmission from the reservoir to the incidental host. Lloyd-Smith et al. (2009) have analysed 442 modelling studies of zoonotic pathogens and concluded that models incorporating spillover transmission are dismayingly rare. More recent models aimed at investigating the dynamics of EIDs by taking into account the spread of the pathogen using multihost processes but disregarding the persistence of the pathogen in the reservoir ( Singh et al., 2013 ) , or by focusing on the dynamics and conditions of persistence of the pathogen between two populations ( Fenton and Pedersen, 2005 ) . Models that have considered an endemic reservoir are disease-specific and do not generate generalizable dynamics ( Chowell et al., 2014; Nieddu et al., 2017 ) . More recently, Singh and Myers (2014) developed a Susceptible-Infected-Recovered (SIR) stochastic model coupled with a constant force of infection. The authors are mostly interested in the effect of population size and its impact on the size of an outbreak. However, this approach does not allow teasing apart the contribution of the incidental host transmission from that of the transmission from the reservoir in modulating the dynamics of zoonotic pathogens.

In this paper, we aim to provide general insights into the dynamics of a zoonotic pathogen (i.e. pathogens classified in stages II-IV) emerging from a reservoir and its ability to propagate in an incidental host. To do so, we developed a continuous time stochas- Fig. 1 . Representation of the classification scheme of pathogens proposed by Wolfe et al. (2007) . A pathogen may evolve from infecting only animals (Stage I) to infecting only humans (Stage V). Each stage corresponds to a specific epidemiological dynamics in the incidental host. Stage II corresponds to few spillovers from animals (e.g. bats) to humans with no possible transmission between humans. Stage III corresponds to few stuttering chains of transmission between humans that go extinct (no outbreaks). Stage IV corresponds to large outbreaks in human population but the pathogen cannot be maintained without the reservoir. tic model that can dissociate the effect of between-host (i.e. direct) transmission from the effect of spillover (i.e. reservoir) transmission. A multi-host process with a reservoir and an incidental host is considered. The epidemiological processes are stochastic, which is particularly relevant in the case of transmission from the reservoir and more realistic because only a small number of individuals are expected to be infected at the beginning of an outbreak. The model makes a number of assumptions. First, the epidemiological processes are much faster than the demographic processes. Second, the pathogen in the reservoir is considered as endemic and might contaminate recurrently the incidental host. Third, an individual cannot become susceptible after having been infected. As a consequence, the total number of susceptible individuals in the incidental host decreases during the epidemic. This is what is expected for an epidemic spreading locally during a short period of time (at the scale of a few thousands individuals during weeks or months, depending on the disease and populations considered). We then harness the model to predict the effects of both spillover transmission and direct transmission on the number and the size of outbreaks. Outbreaks occur when the number of cases of disease increases above the epidemiological threshold. In the case of non emerging infectious diseases, an epidemiological threshold is used to gauge the start of outbreaks. For instance for the seasonal influenza the epidemiological threshold is calculated depending on the incidence of the disease during the previous years ( Tay et al., 2013 ) . In the case of emerging infectious diseases, no incidence is normally expected in the population so from a small number of infected individuals, the outbreak can be considered to spread. We show that, regarding the epidemiological dynamics, the recurrent emergence of the pathogen from the reservoir in the incidental host is as important as the transmission between individuals of the incidental host. We conclude by discussing the implications of these results for the classification of pathogens proposed by Lloyd-Smith et al. (2009) .

A continuous time stochastic Susceptible-Infected-Recovered (SIR) compartmental transmission model ( Kermack and McKendrick, 1927 ) with recurrent introduction of the infection into an incidental host by a reservoir is considered ( Fig. 2 ). Our goal here is not to study a disease in particular but to provide general insights of the reservoir effect on the epidemiological dynamics of the incidental host. The infection is assumed to propagate quickly relatively to other processes such as pathogen evolution and demographic processes. The reservoir is defined as a compartment where the pathogen is persistently maintained, this pathogen is then considered as endemic. The population is fully mixed. An individual can be infected through two types of transmission, from the reservoir by the spillover transmission and by direct contact between individuals. We neglect the possibility for reverse infection from the incidental host to the reservoir.

The incidental host is composed of N individuals. The infection can spillover by contact between the reservoir and the incidental host at rate τ S where S is the number of susceptible individuals and τ is the rate at which an individual becomes infected from the reservoir. In the incidental host, the infection can propagate by direct contact at rate βSI where I is the number of infected individuals and β is the individual rate of infection transmission. An infected individual can recover at rate γ . The propensity of the pathogen to be transmitted between individuals within host is expressed in terms of the basic reproductive ratio of the pathogen, R 0 , which is widely used in epidemiology. R 0 corresponds to the average number of secondary infections produced by an infected individual in an otherwise susceptible population. In a deterministic model, for a pathogen to invade the population, R 0 must be larger than 1 in the absence of reservoir. In a stochastic model, the higher the R 0 the higher the probability for the pathogen to invade the population. In a SIR model, the basic reproductive ratio R 0 equals to βN / γ . Individuals in the recovered compartment do not contribute anymore to the transmission process. Since we assume that demographic processes are slower than epidemic processes, the number of susceptible individuals decreases during the epidemic due to the consumption of susceptible by the infection until the extinction of the susceptible population. In other words, in our model, R 0 will decrease because of the successive spillovers from the reservoir. We expect this to occur especially at short space and time scales (a local population during the course of weeks or months).

To analyse the dynamics in the incidental host, three statistics will be studied (i) the mean number of outbreaks, (ii) the mean size of the recurrent outbreaks during an epidemic and (iii) the 

10 0 0 individuals mean size of the largest outbreak occurring during an epidemic. We consider the appearance of an outbreak when the incidence of the infection exceeds the threshold c and define the maximum size of an outbreak as the largest number of infected individuals during the largest outbreak.

Stochastic simulations. The epidemiological dynamics described previously can be simulated with the following algorithm (simulations were run in C++). The population state is assumed to be known at time t . A total event rate ( ), only depending of the state of the population at time t , is calculated for each iteration.

A) The total event rate of the continuous time stochastic SIR model is given by:

B) The next event time is t = t + δ where δ is exponentially distributed with parameter . C) The next event to occur is randomly chosen: direct transmission, spillover transmission or recover with respective probabilities βSI / , τ S / and γ I / .

We performed stochastic individual-based simulations of the epidemics with spillover transmission, using rates as presented in Fig. 2 . The incidental host is initially ( t = 0 ) composed of 10 0 0 susceptible individuals ( N = S = 10 0 0 ). The infection is considered as endemic in the reservoir. Simulations are stopped when there are no susceptible individuals anymore. An outbreak begins when the number of infected individuals reaches the epidemiological threshold c ( c = 5 infected individuals in the simulations) and ends when there is no infected individuals anymore ( I = 0 ). Stochastic simulations were run for values of the basic reproductive ratio ( R 0 ) ranging from 0 to 10 and of the spillover transmission ( τ ) ranging from 10 −10 to 10 −1 , 10,0 0 0 simulations are performed for each parameter set. All other parameter values are detailed in Table 2 .

Approximation by a branching process. The epidemiological model with recurrent introduction of the infection into an incidental host by a reservoir can be approximated by a branching process with immigration from the reservoir to the incidental host at the beginning of the infectious process (thus assuming that individual "birth" and "death" rates of infected individuals are constant during the starting phase of an outbreak). The individual birth and death rates are respectively βN , the transmission rate and γ , the recovery rate. The immigration rate corresponds to the spillover rate, τ N , at the beginning of the infection. In other words, we assume that the number of susceptibles is N to study the beginning of the infection, which is a good approximation as long as few individuals have been infected. We distinguish between two regimes in the incidental host, the subcritical regime when R 0 < 1 and the supercritical regime when R 0 > 1. We suppose that at time t = 0 a single individual is infected by the spillover transmission.

As illustrated in Fig. 3 , three patterns are observed (i) a stuttering chain of transmission that goes extinct, i.e. infection spreads inefficiently, (corresponding to Stage II in Wolfe's classification, see served are wider for a higher threshold value (compare Fig. 4 a and b), whereas stage III is narrower. When the direct and the spillover transmissions are low, it is more difficult for the infection to reach a higher threshold. Thus, there are more stuttering chains of transmission. In the same way, when the direct or the spillover transmission is high, a large outbreak is observed then some stuttering chains of transmission occur but do not reach the epidemiological threshold.

For both values of epidemiological threshold ( Fig. 4 (a) and (b)), a "bulb" is observed in the Stage III where the direct transmission is high. After the occurrence of a large outbreak, the susceptible population became small. Hence the next excursion is very unlikely to reach the epidemiological threshold. However, a high enough spillover transmission rate is able to counterbalance the small effective R 0 and to produce other outbreaks after the large one. The "bulb" is less pronounced in the case of a higher threshold ( Fig. 4 (b) ) because the susceptible population consumed during the large outbreak is important leading to the failure for the next excursion to reach a high epidemiological threshold.

We aim at approximating the mean number of outbreaks in the case where the spillover transmission rate τ and the reproductive number R 0 are small (subcritical case corresponding to R 0 < 1). The method of approximation is the following: let us denote by S i the number of susceptible individuals at the beginning of the i th excursion. During the i th excursion, we set this number of susceptibles to its initial value S i , and consider that the rate of new infections is βIS i . We thus obtain a branching process with individual birth (infection) rate βS i and individual death (recovery) rate γ .

When there is no more infected individuals, we compute the mean number of recovered individuals produced by this branching process excursion, denoted by E [ K(S i , β, γ )] , and make the approximation that

where E [ K(S i , β, γ )] can be computed and equals (see

(2)

In other words, the initial number of susceptible individuals for the ( i + 1 )th excursion is the initial number of susceptible individuals for the i th excursion minus the mean number of recovered individuals produced during the i th excursion under our branching process approximation. We repeat the procedure for the ( i + 1 )th excursion, and so on, until k satisfies S k > 0 and S k +1 ≤ 0 (no susceptible anymore). In order to be considered as an outbreak, an excursion has to exceed c individuals, where we recall that c is the epidemiological threshold. Under our branching process approximation, the probability for the i th excursion to reach the epidemiological threshold (see Appendix A ) is:

As a consequence, our approximation of the mean number of out-

where S 1 = N, and the S i 's are computed as described in (1) . The mean number of outbreaks computed with the branching process is a good approximation compared to numerical simulations for a small spillover transmission ( 10 −10 ≤ τ ≤ 10 −6 ) ( Fig. 5 ) .

The spillover transmission added in our model introduces the infection recurrently and allows the infection to spread even for a pathogen barely contagious ( R 0 < 1). According to Fig. 5 when R 0 < 1 the number of outbreaks increases when the direct transmission between individuals increases. Indeed, the higher the direct transmission, the higher the probability for the excursions to reach the epidemiological threshold ( c ). The number of outbreaks can be high because when the direct transmission is smaller than 1, the infection spreads inefficiently and does not consume a large number of susceptibles allowing the next excursion to exceed the epidemiological threshold. Fig. 6 shows that the average number of outbreaks is a non-monotonic function of the direct transmission ( R 0 ) and the spillover transmission ( τ ). More precisely, Fig. 6 (b) shows that for intermediate and low values of spillover transmission ( 10 −6 ≤ τ ≤ 10 −4 ), the average number of outbreaks increases until R 0 ∼ 1 then decreases until it reaches one outbreak when the direct transmission is high ( R 0 > 2.5). Moreover, we observed an increasing number of outbreaks with τ when the pathogen is barely contagious.

By contrast, in Fig. 6 (a), we show that the average number of outbreaks decreases when τ becomes large ( τ 5 . 10 −4 ).

The supercritical case ( R 0 > 1) is now considered and the spillover transmission rate ( τ ) is still supposed small. process approximation, respectively. The average number of outbreaks approximated is evaluated when the spillover transmission τ is small. For the numerical simulations, τ = 10 −10 has been chosen. There is a break in the dotted curve (branching process) because our approximation is not valid in the critical regime (when R 0 is close to 1). In this case, two different types of excursions occur in the incidental host: (i) a large outbreak which consumes, with a probability close to one, a large proportion of susceptible individuals and (ii) multiple excursions before and after a large outbreak which each consumes few susceptible individuals. We let O before ( N, β, γ ) and O after ( N, β, γ ) denote the number of outbreaks occurring respectively before and after the large outbreak. Because R 0 > 1, the probability to have one large outbreak is close to one. Hence we make the approximation that one large outbreak occurs during the epidemic, and the total number of outbreaks ( O total 

To be part of outbreaks occurring before the large one, an excursion has to satisfy two conditions (i) to have a size higher than the epidemiological threshold c , and (ii) to be of a size not too large otherwise it would correspond to the large outbreak. To be more precise, this condition will correspond to the fact that the supercritical branching process used to approximate this excursion does not go to infinity. As a consequence, O before ( N, β, γ ) can be approximated by (See Appendix B ) :

To approximate the number of outbreaks after the large outbreak ( O after ( N, β, γ ) ), we need to know how many susceptible individuals remain in the population. The number of susceptibles consumed before the large outbreak is negligible with respect to the number of susceptibles consumed during the large outbreak.

Hence we can consider that the initial state of the large outbreak is N susceptibles, one infected individual and no recovered individual. 

which has one trivial solution ( S = N) and a non-trivial solution with no explicit expression denoted N after ( N, β, γ ) . After the large outbreak, the reproductive ratio for the next excursions, denoted R 0 a f ter , is subcritical ( R 0 a f ter < 1 ) (see Appendix B ) and the number of outbreaks after the large one, denoted O after , can be approximated using Eqs.

(2) -(4) . The branching process approximations of the mean number of outbreaks in the supercritical regime, depicted in Fig. 5 , are close to the mean number of outbreaks computed by numerical simulations when the recurrent infection from the reservoir is small. The number of outbreaks decreases when the pathogen becomes highly contagious to reach one outbreak when R 0 > 2.5. When the infection is introduced in the incidental host by the spillover transmission, the probability to reach the epidemiological threshold depends on the direct transmission between individuals. When the direct transmission increases the infection spreads more efficiently consuming a large number of susceptible individuals allowing few or no other excursion to reach the epidemiological threshold and producing only one outbreak when R 0 > 2.5.

We now focus on the effect of the spillover transmission with a pathogen barely contagious ( R 0 < 1) on the number of outbreaks. We exclude for the sake of simplicity the cases very close to the critical case, that is to say, 1 − R 0 is not too close to 0. Because we consider the subcritical case ( R 0 < 1), the excursions are small and at the beginning of the epidemiological dynamics, we make the approximation that the spillover transmission rate is constant equal to τ N , and the direct transmission rate is equal to βNI . We thus consider a birth and death process with constant immigration rate τ N , individual birth rate βN and individual death rate γ . We are interested in the effect of the parameter τ on the mean number of outbreaks. In particular we aim at estimating the value of τ maximising the mean number of outbreaks, denoted τ opt .

A first quantity which will help giving us an idea of the order of magnitude of the values of τ to be considered is the mean number of infected individuals at large times. This quantity, denoted m I , equals (see for instance Eq. (8.74) in Bailey, 1990 ) :

In particular, when

the mean number of infected individuals is much larger than c , and when on the contrary

the mean number of infected individuals is negligible with respect to c .

Let us first consider the first case ( Eq. (8) ), and choose α > 1 such that

Then we can show (see Appendix C.1 ), that the probability p c that a first infection by the reservoir gives rise to an outbreak (that is to say the number of infected individuals reaches c before 0) is larger than:

Moreover, we can show that if an excursion reaches the level c , it has a probability close to one to lead to a large outbreak consuming a large number of susceptible individuals. Thus only few stuttering chains of transmission will emerge. In Fig. 7 (a) , when τ is large ( τ > 10 −2 thus τ N/ (cγ (1 − R 0 )) ≥ 25 ), only one outbreak is observed because the large number of spillovers prevents the outbreak from dying out. Let us now consider the second case ( Eq. (9) ). Recall that in the case of emerging infectious diseases, the threshold c can be considered as small. Hence we may consider without loss of generality that (9) implies:

In this case, we can prove (see Appendix C.2 ) that the probability that the number of infected individuals is higher than the threshold c is :

Thus the probability for the number of infected individuals to reach the epidemiological threshold c is small under condition (9) . As a consequence, few outbreaks will occur. Indeed, the successive spillovers by the reservoir will produce outbreaks with a small probability, but will nevertheless consume susceptible individuals, until there is no more susceptible in the population. According to Fig. 7 (a) , when a small effect of spillover transmission ( τ < 10 −6 ) and a small reproductive ratio ( R 0 ≤ 0.6) are considered ( τ N/ (cγ (1 − R 0 )) < 2 . 10 −2 ) then the number of outbreaks is small.

In the case of a slightly higher direct transmission rate ( R 0 = 0 . 8 ), each spillover has a non negligible probability to become an outbreak (more precisely 0.12 when c = 5 ) and the number of outbreaks is higher.

We thus predict that the number of outbreaks will tend to be large when the average size of an excursion is close to the epidemiological threshold ( m I c ). These observations allow us to give a first rough upper bound of the optimal value τ opt . Indeed, if the mean number of infected individuals ( E [ I ]) is equal to c , the ra-

represents the variance of the number of infected individuals, see Appendix C.3 ) and whose value belongs to [0.25,1] when c = 5 for the values of R 0 considered, which is large. Moreover in this case the distribution of I is skewed to the right (see Fig. C .10 ). This implies that the number of infected individuals will be larger than c a large fraction of the time producing outbreaks which do not go extinct before a new infection by the reservoir, and thus producing few outbreaks.

Hence we may conclude that τ opt is smaller than the τ leading to a mean number of infected individuals c . For instance for the parameters considered in Fig. 7 (a) , this gives that τ opt is smaller than:

Let us now be more precise on the estimation of τ opt . To this aim, we will apply two results of the theory of branching processes with immigration. The first one, which can be found in Bailey (1990 , Eq. (8.74 )) describes the total infectious lines over the course of the infection, denoted by m :

Notice that m is necessarily larger than 1 as an infection from the reservoir is needed to generate the first infectious line.

The second result is the mean number of infectious lines expected to be present at any time, that is to say in the theory of branching processes the number of distinct immigrants which have descendants alive at a given moment. For large times, this mean number ( M I ) is equal to: Pardoux, 2007 ) , where E [ T 0 ] denotes the mean lifetime of a branching process without immigration, with individual birth rate βN , individual death rate γ , and initial state 1. This expression can be computed explicitly (see Appendix C.3.1 ) and equals:

thus leading to the expression:

We will divide an excursion into m / M I blocks of M I simultaneous infectious lines (thus without immigration). The idea for such an estimation is the following: it is known that if a Poisson process has k jumps during a time interval, the jumps are uniformly distributed during this time interval. As the infections by the reservoir follow approximately a Poisson process with parameter τ N , and we know that in expectation m infections by the reservoir occur before all infected individuals are removed, we divide the epidemic in homogeneous blocks. We choose these blocks to have an initial number of M I infected individuals to allow the use of results on branching processes with immigration. The initial number of infected individuals in each block is thus M I , and as a consequence the infection has a probability

to reach the threshold c (see Appendix C.5 ). Hence the probability for the whole excursion to reach the threshold c can be approximated by

We want this probability to be not too close to 0, otherwise most susceptible individuals would be consumed without giving rise to an outbreak. We also want this probability to be not too close to 1. Indeed, as we have shown in the beginning of this section, this would correspond to a case where τ N / γ is much larger than c(1 − R 0 ) and once the infected number of individuals has reached the value c it would be very likely to reach a large value and consume a large number of susceptible individuals. As a consequence, we would have at the limit only one large outbreak. We thus choose to equalize this probability to one half to get an estimation of τ opt .

Notice that this choice is arbitrary but has only a small effect on the final results. For instance, a choice of 0.3 or 0.7 would give very close results. The most important is to stay away from 0 and 1. As a conclusion, τ opt is estimated as the unique solution to:

The unicity of the solution is proved in Appendix C.5 . Fig. 7 (b) presents the values of τ maximising the number of outbreaks and their estimations (dots) obtained by the branching process approximations. The estimates derived under the branching process approximation give good results, with error ranging from 3 to 33 % regardless the value of the epidemiological threshold ( Table 3 ) .

Values of the optimal spillover transmission from numerical simulations and estimations. The optimal spillover transmission is calculated for R 0 being equal to 0.2, 0.4, 0.6, and 0.8. We present the values of the optimal spillover for two values of epidemiological threshold c = 5 (bold) and c = 10 (not bold). The errors are ranging from 3 to 33% with a mean error of 15%. To get the estimation of τ opt we have made several approximations. First we have considered that the spillover rate by the reservoir is constant equal to τ N , whereas it is decreasing and equals to τ S , and that the rate of direct transmission due to infected individuals in the population is βNI and not βSI . We believe that these approximations are reasonable because the probability for an excursion to reach the threshold decreases with the consumption of susceptible individuals, and as a consequence, most of the outbreaks will occur at the beginning of the process. However, the real τ opt should be a little bit higher than the one we estimate, to counterbalance the fact that the real infection rates (by the reservoir and the infected individuals) are smaller than the one we use in our calculations. This may explain why in most of the cases we underestimate the real τ opt (see Table 3 ).

The following approximation we made is the decomposition of the excursions into blocks with an initial number of individuals M I . In the real process there are no simultaneous infections by the reservoir. However this approximation allows to take into account previous infections by the reservoir whose infectious lines are still present.

When the ratio, V ar[ I ] / E 2 [ I ] = τ N/γ , is small (see Appendix C.3 ), the value of I stays close to its expectation and few outbreaks occur, as the number of infected individuals rarely reaches 0. For instance, for R 0 = 0 . 2 and τ = 4 . 9 . 10 −4 , τ N / γ ∼ 0.2. As decreasing τ increases this ratio, this could explain why we overestimate τ opt for small values of R 0 (because smaller values of R 0 necessitate higher values of τ to get the same probability to reach the threshold c ). 

During the epidemic, a large outbreak can occur depending on the value of the direct transmission ( R 0 ) and the spillover transmission ( τ ) and corresponds to the largest number of infected individuals. To analyse the effect of the recurrent emergence of the pathogen on the size of the largest outbreak, we model the largest outbreak by a SIR deterministic model with a spillover transmission:

Since no explicit expression of the size of the outbreak can be obtained with the deterministic model, we estimated it using numerical analyses. Fig. 8 shows that the maximal number of infected individuals during the largest outbreak increases with the direct transmission ( R 0 ) and the spillover transmission ( τ ). When the direct transmission ( R 0 ) is small, the size of the largest outbreak can differ by orders of magnitude with varying spillover transmission ( τ ). Furthermore, a large outbreak can be observed for a pathogen barely contagious ( R 0 < 1) when the recurrent emergence of the pathogen is high ( τ 10 −3 ).

Zoonotic pathogens constitute one of the most pressing concerns with regards to future emerging diseases, but studies investigating the importance of the role of animal reservoirs for the epidemiological dynamics of infectious diseases are lacking. Indeed, most theoretical works only consider pathogen transmission between conspecifics for predicting disease epidemiology. Here, we build a continuous time stochastic SIR model to consider the statistical process underlying a spillover transmission, i.e. transmission from an animal reservoir to a host. We analyse the model to predict the number and the size of outbreaks as a function of both the spillover transmission and within host. The model shows that spillover transmission influences the epidemiological dynamics as much as the transmission by direct contact between individuals. Three different dynamics are observed, ranging from the absence of outbreaks to a single large outbreak. The findings have implications for (1) modelling the dynamics of EIDs, (2) understanding the occurrence of outbreaks in the case of pathogens that are barely contagious and (3) control strategies.

In our results, the appearance of outbreaks depends on both the transmission from the reservoir and the direct transmission between individuals. Generally, the occurrence of epidemics in humans is attributed to the ability of the pathogen to propagate between individuals. In the case of a single-host process, the notion of the basic reproductive ratio R 0 seems sufficient to evaluate the spread of the pathogen in a population entirely composed of susceptible individuals. In EIDs, R 0 is also used to gauge the risk of pandemics. In this way, Lloyd-Smith et al. (2009) delineate the three stages identified for a zoonotic pathogen ( Wolfe et al., 2007 ) by using the ability of the pathogen to spread between individuals. Each stage corresponds to a specific epidemiological dynamics ranging from a non-contagious pathogen making an outbreak impossible (Stage II, R 0 = 0 ) to a barely contagious pathogen with few outbreaks and stuttering chains of transmission (Stage III, R 0 < 1) to a contagious pathogen making a large outbreak possible (Stage IV, R 0 > 1). The aim of the Wolfe's classification is to establish each stage at which a zoonotic pathogen may evolve to be adapted to human transmission only, in order to identify pathogens at potential risk of pandemics. However, in our model, by taking into account the recurrent emergence of the pathogen from the reservoir, the three dynamics that define the three stages will depend on both the spillover transmission and the direct transmission of the pathogen between individuals. The results suggest that in the case of pathogen spilling recurrently over an incidental host, the direct transmission should not be the only parameter to consider.

The presence of a reservoir and its associated recurrent spillovers dramatically impact the epidemiological dynamics of infectious diseases in the incidental host. Without transmission from the reservoir, the probability to have an outbreak when the pathogen is barely contagious only depends on the direct transmission between individuals, and the outbreak rapidly goes extinct. By contrast, the results show that the recurrent emergence of the pathogen from a reservoir increases the probability to observe an outbreak. Spillover transmission enhances the probability to both observe longer chains of transmission and reach the epidemiological threshold (i.e. threshold from which an outbreak is considered) even for a pathogen barely contagious. Moreover, this coupling model (reservoir-human transmission) allows the appearance of multiple outbreaks depending on both the ability of the pathogen to propagate in the population and the transmission from the reservoir. Zoonotic pathogens such as MERS, Ebola or Nipah are poorly transmitted between individuals ( R 0 estimated to be less than 1) ( Althaus, 2014; Chowell et al., 2014; Luby et al., 2009; Zumla et al., 2015 ) yet outbreaks of dozens/hundreds/thousands of infected individuals are observed. We argue that, as suggested by our model, the human epidemic caused by EIDs could be due to a recurrent spillover from an animal reservoir.

In the case of zoonotic pathogens, it is of primary importance to distinguish between primary cases (i.e. individuals infected from the reservoir) and secondary cases (i.e. individuals infected from another infected individual) to specify which control strategies to implement and how to optimize public health resources. According to the stochastic SIR model coupled with a reservoir analysed here, the same dynamics can be observed depending on the relative contribution of the transmission from the reservoir and the direct transmission by contact with an infected individual (see Fig. 4 ). For example, a large outbreak is observed either for a high spillover transmission or for a high direct transmission. The proposed stochastic model makes it possible to understand the effects of the infection from the reservoir or from direct transmission on the epidemiological dynamics in an incidental host when empirically this distinction is difficult. Empirically, the origin of the infection is established by determining the contact patterns of in-fected individuals during the incubation period. Thereafter, the role of control programs implemented could be evaluated in order to determine better strategies.

We have considered that the reservoir is a unique population in which the pathogen can persist, which is a simplifying assumption. The pathogen is then endemic in the reservoir and the reservoir has a constant force of infection on the incidental host. The reservoir can be seen as an ecological system comprising several species or populations in order to maintain the pathogen indefinitely ( Haydon et al., 2002 ) . For example, bat and dromedary camel ( Camelus Dromedarius ) populations are involved in the persistence of MERS-CoV and in the transmission to human populations ( Sabir et al., 2016 ) . In these cases, the assumptions of a constant force of infection can be valid because the pathogen is endemic. However, zoonotic pathogens can spill over multiple incidental hosts and they can infect each other. In the case of the Ebola virus, which infects multiple incidental hosts such as apes, gorillas and monkeys ( Ghazanfar et al., 2015 ) , the principal mode of contamination of the human population is the transmission from non-human primate populations. Moreover, the contact patterns between animals and humans is one of the most important risk factors in the emergence of avian influenza outbreaks ( Meyer et al., 2017 ) . These different epidemiological dynamics with transmission either from the reservoir or from other incidental hosts can largely impact the dynamics of infection observed in the human population, and the investigations of those effects can enhance our understanding of zoonotic pathogens dynamics.

In our model, we make a second simplifying assumption by considering that the infection propagates quickly relatively to other processes such as pathogen evolution and demographic processes. This assumption can be not valid in the case of low emergence of the pathogen from the reservoir. Indeed, the time between two spillovers can be long and makes the evolution of the pathogen possible inside the reservoir. Moreover, during the time between two spillovers, the demography in the incidental host can vary and impact the propagation of the pathogen. In the case of low spillover transmission in the incidental host, the effect of both pathogen evolution and demographic processes can be a topic for future research on the epidemiological dynamics of emerging infectious diseases.

In this paper, we have argued that the conventional way for modelling the epidemiological dynamics of endemic pathogens in an incidental host should be enhanced to account for spillover transmission in addition to conspecifics transmission. We have shown that our continuous time stochastic SIR model with a reservoir produces similar dynamics to those found empirically (see the classification scheme for pathogens from Wolfe et al., 2007 ) . This model can be used to better understand the ways in which EIDs transmission routes impact disease dynamics.

In this appendix, we derive results on the branching process approximation stated in Section 3 . The main idea of this approximation is the following: when the epidemiological process is subcritical ( R 0 < 1), an excursion will modify the state of a small number of individuals with respect to the total population size. During the i th excursions, the direct transmission rate βSI will stay close to βS i I where S i denotes the number of susceptibles at the beginning of the i th excursion. Hence, if we are interested in the infected population, the rate βS i I can be seen as a constant individual birth rate βS i . Similarly, γ I , which is the rate at which an individual in the population recovers, can be interpreted as a constant individual death rate γ in the population of infected individuals.

In this section, we focus on the number of outbreaks when R 0 < 1 and when the rate of introduction of the infection by the reservoir is small ( τ 1). That is to say, we consider that each introduction of the infection by the reservoir occurs after the end of the previous excursion created by the previous introduction of the infection by the reservoir. According to Eq. (12) , this approximation is valid as long as the ratio τ / β is small. We first approximate the mean number of susceptible individuals consumed by an excursion. Let us consider a subcritical branching process with individual birth rate βN and individual death rate γ . As this process is subcritical, we know that the excursion will die out in a finite time and produce a finite number of individuals. Then from Britton and Pardoux (2018) or Van Der Hofstad (2016) , if we denote by K [ N, β, γ ] the total number of individuals born during the lifetime of this branching process (counting the initial individual), we know that:

where P denotes a probability, and hence

where E is the expectation.

By definition, an excursion is considered as an outbreak only if the maximal number of individuals infected at the same time during this excursion is larger than an epidemiological threshold that we have denoted by c . Hence in order to approximate the number of outbreaks we still have to compute the probability for an excursion to be an outbreak. This is a classical result in branching process theory which can be found in Athreya and Ney (1972) for instance. P (N, β, γ ) : = P ( more than c individuals infected at a given time ) = γ /βN − 1 (γ /βN) c − 1 .

(A.1)

With these results in hands, the method to approximate the mean number of outbreaks is the following: the probability that the first excursion is an outbreak is Fig. A.1 . Probability for an excursion to be of size k . The situation is indicated for a reproductive ratio ( R 0 ) varying from 0.2 to 0.8. The number of susceptibles at the beginning of the second excursion is approximated by

The second excursion has a probability γ /βS 2 − 1 (γ /βS 2 ) c − 1 to be an outbreak. The number of susceptibles at the beginning of the third excursion is approximated by

and the third excursion has a probability γ /βS 3 − 1 (γ /βS 3 ) c − 1 to be an outbreak. The procedure is iterated as long as there is still a positive number of susceptible individuals. This gives 4 .

We now focus on the case R 0 = βN/γ > 1 . In this case the approximating branching process is supercritical and go to infinity with a positive probability. In the case where the epidemic process describes small excursions, the branching process approximation is still valid, but in the case where it describes a large excursion, then a large fraction of susceptible individuals is consumed and the branching approximation is not valid anymore. However, as all the quantities (susceptible, infected and recovered individuals) are large, a mean field approximation is a good approximation of the process. Here the mean field approximation is the deterministic SIR process, whose dynamics is given by:

Let us first focus on the small excursions occurring before the large one.

As they are small, they can be approximated by a branching process. Here, unlike in the previous section, the approximating branching process Z is supercritical, as βN > γ . We compute its probability to drift to infinity:

As we will see, a supercritical branching process with individual birth rate βN and individual death rate γ conditioned to go extinct has the same law as a subcritical branching process with individual birth rate γ and individual death rate βN . Indeed, if we denote by Z n the successive values of this branching process, we get for every couple of natural numbers ( n, k ):

where P (A | B ) denotes the probability of the event A when B is realised. We used again in this series of equalities classical results on branching processes that can be found in Athreya and Ney (1972) .

As a consequence, if we denote by G [ N, β, γ ] the number of susceptible individuals consumed by the excursion of a supercritical branching process with individual birth rate βN and individual death rate γ conditioned to go extinct, we get:

, and the probability for this excursion to have a size bigger than the epidemiological threshold c is

As the number of susceptible individuals stays large until the large excursion occurs, we may keep N as the initial number of susceptibles at the beginning of the excursions instead of replacing it by their mean value, as we have done in the previous section.

The different quantities we have just computed allow us to approximate the number of small excursions before the large excursion: in expectation, we have

As k and l ( k ) are small with respect to N , this can be approximated by

Finally, notice that in (B.1) , S is a decreasing quantity, and I is a non-negative quantity, which varies continuously. Hence ˙ I = I(βS − γ ) has to be negative before I hits 0. As a consequence,

This ensures that the epidemic is subcritical after the large outbreak.

In this section, we focus on the effect of the reservoir transmission rate (parameter τ ) on the number of outbreaks when the infection is subcritical ( R 0 < 1). The idea is the following: first, as excursions of subcritical branching processes are small, we can make the approximation that, at the beginning, the infection rate by the reservoir is constant equal to τ N , and that the direct transmission rate is equal to βNI . Making this approximation allows us to handle the two processes of infection (by contact and by the reservoir) independently, and to use known results on branching processes with immigration.

Recall that I denotes the number of infected individuals. We first assume (8) and prove inequality (10) . Let us choose α > 1 such

Then for any 1 ≤ k ≤ c − 1 , the jump rates of the process I are:

This implies that once one individual is infected, the probability for the number of simultaneously infected individuals to reach c before the recovery of all infected individuals is larger than the probability that a birth and death process with initial state 1 and birth ( b ) and death ( d ) rates satisfying b d

reaches the state c . Applying (A.1) , we deduce that this probability p c C.3 satisfies: 

. This proves (10) .

Moreover, as α has been taken large, the infectious process stays supercritical (in the sense that the next event is more likely to be an infection than a recovery) until a size k satisfying

and thus if the number of infected individuals reaches c it is likely to reach a large value and consume a large number of susceptible individuals.

As we approximate the infection process by a branching process with constant immigration, the law of I under this approximation converges to a well-known law, provided by Eq. (8.75) in Bailey (1990) :

From this law, we can deduce the probability for I to be equal to any integer k :

and thus, the probability for I to be larger than c is

Recall that we assumed that

We thus get the following inequality

Now let us notice that for any i ≥ 1:

1 − R 0 2

We thus get

where to get the last inequality we computed the maximum of the function x → xa x −1 for a = (1 + R 0 ) / 2 . As a conclusion, for a fixed R 0 < 1 and a small enough τ , the probability for the number of infected individuals at a given time to be higher than c is bounded by a function of R 0 time τ N / (cγ (1 − R 0 )) . This probability is thus small when the last term is small.

Recall Eq. (C.1) . It allows us to compute the variance of I , as follows:

In particular,

E 2 [ I] = τ N γ .

In this section we provide an expression for the term E [ T 0 ] which appears in the definition of M I (see Eq. (13) ). This expression derives from the following equality, which can be found in Athreya and Ney (1972) . Let t ≥ 0 and T 0 denotes the extinction time of the excursion of a branching process, with one individual at time 0. Then we have

This allows one to compute the expectation of T 0 as follows:

where we made an integration by parts. Thus

Recall that if we consider a branching process with individual birth rate βN , individual death rate γ , and initial state k ≤ c , the probability for this process to reach the size c is Athreya and Ney, 1972 for instance). We use this result to approximate the probability for a block of the excursion to reach the threshold c by

where we recall that M I is the mean number of simultaneous excursions generated by different infections from the reservoir. Notice that this is an approximation, as M I is not necessarily an integer.

We end this appendix with the proof of the unicity of the solution to (14) . To simplify the notations, we introduce the function F , which at x associates:

First we notice that F is only defined for x ≤ c . Otherwise the term in brackets would be negative. Second, notice that if x ≤ 1 and for any c ≥ 2 and R 0 < 1,

This shows that if x ≤ 1, F ( x ) > 1/2. We now determine the sign of F ( x ) for x belonging to the interval [1, c ]. A direct computation gives:

As the two logarithms are negative for x belonging to (1, c ), we deduce that F ( x ) < 0 for x belonging to (1, c ). As F (1) > 1/2 and F (c) = 0 , we conclude that F (x ) = 1 / 2 has a unique solution on the interval (1, c ). This is equivalent to the fact that (C.2) has a unique solution τ opt which belongs to

This ends the proof.

Estimating the reproduction number of ebola virus (ebov) during the 2014 outbreak in west africa

What it takes to be a reservoir

When is a reservoir not a reservoir?

Branching Processes

The Elements of Stochastic Processes with Applications to the Natural Sciences

Centre for Disease Control and Prevention

Synthesizing data and models for the spread of mers-cov, 2013: key role of index cases and hospital transmission

Ebolavirus evolution: past and present

Community epidemiology framework for classifying disease threats

Ebola, the killer virus

Identifying reservoirs of infection: a conceptual and practical challenge

Ebola between outbreaks: intensied ebola hemorrhagic fever surveillance in the democratic republic of the congo

Global trends in emerging infectious diseases

Impacts of biodiversity on the emergence and transmission of infectious diseases

A contribution to the mathematical theory of epidemics

Epidemic dynamics at the human-animal interface

Recurrent zoonotic transmission of nipah virus into humans

Movement and contact patterns of long-distance free-grazing ducks and avian influenza persistence in vietnam

Emerging infectious diseases: threats to human health and global stability

Human ecology in pathogenic landscapes: two hypotheses on how land use change drives viral emergence

Extinction pathways and outbreak vulnerability in a stochastic ebola model

Processus de markov et applications: algorithmes, réseaux, génome et finance: cours et exercices corrigés

Co-circulation of three camel coronavirus species and recombination of mers-covs in saudi arabia

Outbreak statistics and scaling laws for externally driven epidemics

The structure of infectious disease outbreaks across the animal-human interface

Exploring a proposed who method to determine thresholds for seasonal influenza surveillance

Risk factors for human disease emergence

Random Graphs and Complex Networks

Origins of major human infectious diseases

Biological features of novel avian influenza a (h7n9) virus

Middle east respiratory syndrome

The authors have been supported by the "Chair Modélisation Mathématique et Biodiversité" of Veolia Environnement-Ecole Polytechnique-Museum National d'Histoire Naturelle-Fondation X, France.

outbreaks. Now we focus on the large excursion. We use Eq. (B.1) to approximate its dynamics. This equation is well-known, and it is easy to obtain the equation satisfied by the final number of susceptible individuals: from (B.1)In particular, if we are interested in the time T f when there is no more infected individual and we suppose that at time 0 there is only one infected individual we getThat is to say, S ( T f ) and S (0) are related by the equationRigorously, the value of S (0) depends on the number of susceptible individuals consumed by the small excursions before the large excursion, but we have seen that this number is small compared to the population size, N . Hence the number of susceptible individuals remaining after the large excursion can be approximated by the smallest solution of:Notice that it is easy to have an idea of the error caused by a small variation of the initial state. Indeed, if we denote by S f the smallest solution of (B.2) and by S f − l(k ) the solution when
