Stochastic Epidemics in Growing Populations

Consider a uniformly mixing population which grows as a super-critical linear birth and death process. At some time an infectious disease (of SIR or SEIR type) is introduced by one individual being infected from outside. It is shown that three different scenarios may occur: (i) an epidemic never takes off, (ii) an epidemic gets going and grows but at a slower rate than the community thus still being negligible in terms of population fractions, or (iii) an epidemic takes off and grows quicker than the community eventually leading to an endemic equilibrium. Depending on the parameter values, either scenario (i) is the only possibility, both scenarios (i) and (ii) are possible, or scenarios (i) and (iii) are possible.


Introduction
Consider a population in which an infectious disease is introduced. The question of what might happen is the topic of many scientific papers, relying on mathematical or statistical models, or simulations, and often also fitted to specific populations and/or diseases.
In the current paper, we use stochastic models to answer this type of questions for a dynamic population model. We explicitly include that also the population changes randomly over time: old individuals die and new are born, in such a way that the population, on average, grows. Our population is modeled by a super-critical linear birth and death process. An infectious disease is then introduced into this population. The disease is of SIR (Susceptible-Infectious-Removed) type (e.g., Diekmann et al. 2012) meaning that individuals are at first Susceptible. If they get infected they become Infectious, and after some time in the infectious stage they Recover and become immune and stay so forever. In Sect. 4, we extend the result to an SEIR epidemic where individuals upon infection at first become latent (Exposed but not yet infectious) for some time before they enter the Infectious state.
We study the simplified model where individuals mix (and hence also infect) uniformly in the community. We prove that there are three qualitatively different scenarios that may happen. Either the disease outbreak dies out quickly. If not, there may be a big outbreak infecting more and more people, but still infecting only a negligible minority of the (growing) population. The third scenario is where the epidemic outbreak grows faster than the population, leading to an endemic situation in which the fraction of infectives fluctuates around some fixed level, the endemic level. Depending on the parameter values either only the first scenario is possible, the first and second, or the first and third scenarios.
Infectious diseases in growing populations have been modeled by several authors. Thieme (1992) and Li et al. (1999) study fatal deterministic epidemic models allowing the population to be growing. Diekmann et al. (2012, Sect. 4.4) study a model very similar to the current SIR model, also allowing for disease fatalities, the difference being that we here look at a stochastic model for a finite population. Breban et al. (2005) consider approximations of epidemics taking place on networks in growing populations.
The rest of the paper is structured as follows. In Sect. 2, the growing population model and the SIR epidemic model are defined. In Sect. 3, we state the results for the model and give the proofs. In Sect. 4, we discuss how to extend the results to an SEIR epidemic model, and discuss interesting future problems.

The Dynamic Population Model
The population model is defined to be a super-critical Markovian branching process, or more specifically, a super-critical linear birth-and-death (B-D) process. More precisely, the life lengths of individuals are exponentially distributed with death rate μ, and during the life each individual gives birth to new individuals according to a Poisson process with intensity λ, all random quantities mentioned being independent of each other. We assume that λ > μ, i.e., that the process is super-critical. An alternative description of the process, the B-D description, is by the rates at which the process increases by one and decreases by one. If N(t) denotes the number of individuals at t, then the rate of increase is λN (t) and the rate of decrease is μN (t). The fact that N(t) appears linearly in both rates makes it a linear B-D process, and the fact that the net increase/decrease is λ − μ > 0 (by assumption) makes the process super-critical. The starting state is N(0) = n for some strictly positive integer n.
It is well-known (e.g., Haccou et al. 2005;Jagers 1975) that N(t) has positive probability of growing to infinity, and if it does, it grows at the rate N(t) ∼ e (λ−μ)t . By this we mean that if n = 1 then N(t)/e (λ−μ)t converges to a random variable X with an atom of probability μ/λ at X = 0, and conditioned on not being 0, X is exponentially distributed with intensity parameter 1 − μ/λ. In what follows, we are only interested in the part of the sample space where this happens. For that reason we assume N(0) = n is large (whenever n is important, notations are equipped with an n-index). For any fixed ε we can, in fact, choose n 0 large enough such that for any n ≥ n 0 we have Further, it is also known (Jagers 1975, Corollary 6.8.2) that the age distribution among living individuals converges to an exponential distributed random variable with mean 1/λ. More precisely, the age A of a randomly selected individual converges (as t → ∞) to an exponential distributed random variable with mean 1/λ, so the death rate does not affect the asymptotic age distribution.

The Uniform Markovian Dynamic SIR Epidemic Model
We now define a uniformly mixing Markovian epidemic model "on" the dynamic population model described in the previous section. We assume an SIR (Susceptible-Infectious-Removed) epidemic and extend this to an SEIR epidemic model, also having an Exposed/latent state when infected and prior to the infectious state, in Sect. 4, cf. Diekmann et al. (2012).
Initially (t = 0) the population consists of N(0) = n individuals. At this time, an SIR infectious disease is introduced by infecting one uniformly chosen individual. The epidemic is Markovian: An infectious individual remains infectious for an exponential time with rate δ (unless, of course, it dies before). When the infectious period stops, he/she recovers and becomes immune for the rest of his/her life. During the infectious period, the individual has infectious contacts randomly in time according to a homogeneous Poisson process with rate κ. (Note that κ is assumed to be independent of n, so when increasing n individuals do not meet others at higher rates. If κ were to grow with n, the effect would be that R 0 defined below would tend to infinity.) Each time this happens the contacted person is randomly selected among all living individuals. If the contacted person is susceptible, he/she becomes infected and infectious-otherwise the contact has no effect. There is no vertical transmission (i.e., all new born individuals are susceptible) and the disease has no effect on the death rate which hence still equals μ.
Let Z(t) = (S(t), I (t), R(t)) respectively denote the number of susceptible, infectious and recovered individuals at time t. The epidemic is initiated at Z(0) = (S(0), I (0), R(0)) = (n − 1, 1, 0) (the index n is added when the initial population size is important). The possible events and their rates, when currently in state Z(t) = (S(t), I (t), R(t)) = (s, i, r) = z, are given in Table 1.
To summarize, the population model has two parameters, the birth rate λ and the death rate μ (where λ > μ is assumed), and the epidemic has two parameters: the infection rate κ and the recovery rate δ.

Results
We divide the results into the initial phase of an epidemic, where the fraction infected is small, and those where the epidemic grows to an endemic state.

Initial Phase of the Epidemic
We now study the epidemic during the initial phase (when the fraction infected is still small). If we focus on the number of infectives I (t), we see from the table above that I (t) increases at rate κI (t) · S(t)/N(t) (due to infection) and decreases at rate (δ + μ)I (t) (due to recovery or death). As a consequence, when the fraction infected is still small, or equivalently, the fraction susceptible S(t)/N(t) is close to 1, then the number of infectives I (t) behaves almost like a linear branching process with birth rate κI (t) and death rate (δ + μ)I (t). In what follows, we will make this approximation precise. Since S(t)/N(t) ≤ 1, we can dominate the number of infectives by an "upper" linear branching process I (U ) (t) with I (U ) (0) = 1, and having linear birth rate κI (U ) (t) and linear death rate (δ + μ)I (U ) (t). Since the birth rate for I (U ) (t) is always as large as that of I (t) and the death rates are the same, the two processes can be coupled such that I (U ) (t) ≥ I (t) almost surely for all t using methods very similar to those of Ball and Donnelly (1995).
Obtaining a "good" lower bound is somewhat more involved. We do this for the process I (t) up until the first time that S(t)/N(t) ≤ 1 − ε 0 , where ε 0 > 0 can be chosen arbitrarily small. This method is inspired by Durrett (2007), see also Whittle (1955). Before this time-point I (t) can be bounded from below by the linear B-D process I (L) (t) having the same death rate (δ + μ)I (L) (t), but with birth rate κ(1 − ε 0 )I (L) (t). All three processes can be defined on the same probability space (using the same random life-lengths) such that I (L) The number of infectives in the epidemic process is, at least until a fraction ε 0 of all living individuals have been infected, sandwiched between the two linear B-D processes, both having individual death rate μ + δ, and individual birth rate κ(1 − ε 0 ) and κ, respectively. We use this to prove what will happen during the early stages of the epidemic in the theorem below. First, we define two important quantities for the epidemic/branching processes. The first is the basic reproduction number, which is defined as the expected number of infectious contacts an infected individual has. For the current model this number equals This follows since an individual has contacts at rate κ during the infectious period which lasts for a mean duration of 1/(δ + μ); the individual stops being infectious due to recovery or death. The Malthusian parameter α is defined as the exponential growth/decay rate the epidemic/branching process has. It is given as the solution to where c(t) is the expected rate at which an individual has infectious contacts t time units after it was infected (Jagers 1975). For the current model, the contact rate is κ and the contact is infectious if the individual is still infected, which it is with probability e −(δ+μ)t . Putting this into the equation above and solving for α gives the Malthusian parameter being equal to From Eqs. (1) and (2) it follows that the basic reproduction number R 0 exceeds the value of 1 if and only if α exceeds 0. The interpretation of this general rule is that a major outbreak can occur if and only if R 0 > 1 and this also implies that the epidemic will in expectation grow exponentially. We now state our main theorem concerning the early stage of an outbreak. What may happen depends on the sign of α = κ − (δ + μ) and on its relation to the growth rate λ − μ of the population.
Remark What happens after the initial growth in case (iv) is treated in the next section.
Proof (i) The first statement follows from the fact shown above the theorem that I n (t) can be sandwiched between two linear birth and death processes having arbitrarily close behavior on every finite time interval.
(ii) It follows from standard birth-and-death process theory that I (U ) is sub-critical and hence goes extinct with probability 1. The coupling I n (t) ≤ I (U ) (t) applies forever, which then gives the result.
(iii) Since we may consider arbitrarily large n and hence arbitrarily small ε 0 , both bounding B-D-processes are super-critical. Their probabilities of extinction (i.e., minor outbreak) are obtained by first deriving the distribution of the number of infectious contacts that one infected has, i.e., the offspring distribution, and then deriving the extinction probability of the corresponding branching process using theory for branching processes (Jagers 1975). In the current situation, we have exponential duration of the infectious period, the infectious period stops at rate δ + μ, and during this period an infective infects others at rate κ and κ(1 − ε), respectively. As is wellknown (Diekmann et al. 2012), the offspring distribution is then geometric, and the minor outbreak probability equals (δ + μ)/κ and (δ + μ)/(κ(1 − ε 0 )), respectively, but, since ε 0 may be made arbitrarily small, the limiting probability of a minor outbreak equals (δ + μ)/κ. With the remaining probability the two B-D-processes grow exponentially at rate κ − (δ + μ) and κ(1 − ε 0 ) − (δ + μ), respectively. Since the population grows at a higher exponential rate λ − μ and we start with arbitrarily large population N n (0) = n, the ratio S n (t)/N n (t) = 1 − (I n (t) + R n (t))/N n (t) will with arbitrarily large probability never go below 1 − ε 0 (so the coupling never breaks down) and eventually I n (t)/N n (t) → 0.
(iv) The initial statements are identical to those in (iii). The only difference is that the population growth rate of I n (t) now is greater than that of N n (t) (at least if κ > δ + λ), so eventually S n (t)/N n (t) will be below 1 − ε 0 , and the coupling breaks down. In the next section, we show that in that case I n (t) → ∞ as t → ∞.
Case (iii) of the theorem might bring up the following question: Can we couple I (U ) (t) and I n (t) such that I (U ) (t) − I n (t) = 0 hold forever with positive probability? Using the theory as presented in Ball and Donnelly (1995), we note that I (U ) (t) − I n (t) = 0 forever with probability t∈T S n (t−)/N n (t−), where T is the set of times at which I (U )

(t) increases by one. This probability is positive if and only if
From the theory of branching processes and part (iii) of the theorem, we know that if the epidemic does not go extinct, then e −(λ−μ)t N n (t), e −αt I n (t) and e −αt (I n (t) + R n (t)) converges with probability 1, to positive random variables (say W 1 , W 2 and W 3 , respectively). Heuristic arguments now give that for every n and a large enough constant c, t∈T (I n (t−)+R n (t−))/N n (t−) ≤ c ∞ 0 (κW 2 e αt W 3 e αt /W 3 e (λ−μ)t ) dt. This sum is finite if α 2 < λ − μ, which suggests that the event I (U ) (t) − I n (t) = 0 forever has positive probability only on the part of the parameter domain satisfying 0 < α < λ − μ.

Endemicity
In this section, we still assume that λ − μ > 0 (super-critical population process) and we only consider the case α > λ − μ (the case α = λ − μ is omitted), so that the epidemic, in case it takes off, grows at a higher rate than the population. We focus only on the situation where I n (t) and I (U ) (t) tend to infinity; the situation where they go extinct being treated in the previous section.
As mentioned in Sect. 2.1, the population size N n (t) will be close to the deterministic function n(t) = ne (λ−μ)t in the sense that if n = N n (0) is chosen large enough, then |N n (t)/n(t) − 1| < ε for all t with probability arbitrarily close to 1 (see Sect. 2). Since the population grows we introduce the "proportion"-processesS n (t

) = S(t)/n(t),Ī n (t) = I (t)/n(t) andR n (t) = R(t)/n(t), andZ n (t) = (S n (t),Ī n (t),R n (t)).
Below we show that the vector-valued stochastic processZ n (t) converges to a deterministic vector-processz(t). In order to obtain a non-trivial deterministic process, it is necessary to start the proportion process with a positive, albeit small, fraction infectives. For this reason we assume thatZ n (0) = (S n (0),Ī n (0),R n (0)) = n −1 (S(0), I (0), R(0)) = (1 − ε 1 − ε 2 , ε 1 , ε 2 ) for some small but otherwise arbitrary ε 1 , ε 2 > 0. This starting point is arbitrarily chosen and we do not claim that it is exactly here that the stochastic epidemic process will pass through when the number of infectives grows in comparison with the population. In what follows, we show that the starting point has hardly any effect on the state of the process after a long time as long asS n (0) is close to 1 andĪ n (0) andR n (0) are close to 0.
The proof thatZ n (·) converges toz(·) uses methods from Ethier and Kurtz (2005, Chap. 5), see also Andersson and Britton (2000, Chap. 5) applying the theory to epidemic processes, with the difference that now a time-inhomogeneous normalization is used. The epidemic process is defined in terms of (stochastic) rates (see Table 1). As in Ethier and Kurtz (2005) this implies that the epidemic process can equivalently be defined using Poisson processes. Let, for each of the six types of jumps (specified by ) defined in Table 1, Y (·) denote independent standard (intensity 1) Poisson processes, and letŶ (·) be the corresponding centered processes, soŶ (t) = Y (t) − t. We start by writing the original epidemic process Z(t) = (S(t), I (t), R(t)) using Poisson processes: This follows from Ethier and Kurtz (2005) but can also be explained heuristically. The process starts in the correct point, and at any given time point t the rate of the different jumps are as defined in Table 1 and whenever a particular jump occurs the process Z(t) = (S(t), I (t), R(t)) changes with the vector . The different Poisson processes were defined independently, the dependencies lie in for how long they are observed. We now rewrite Eq. (3) in terms ofZ n (t). First, we observe that β (Z(u)) = n(u)β (Z n (u)) for all , except = (−1, +1, 0) and = (+1, 0, 0) for which we have β (−1,+1,0) (Z(u)) = (n 2 (u)/N (u))β (−1,+1,0) (Z n (u)) and β (1,0,0) (Z(u)) = N(u)β (1,0,0) (Z n (u)). However, since N(u)/n(u) is arbitrarily close to 1 with arbitrarily large probability, the same relation approximately holds also for these (the error term is neglected in what follows). Using this we obtain The middle term on the right will become small as n grows (recall that n(t) = ne (λ−μ)t ) because a centered Poisson process divided by something proportional to the mean of the non-centered process becomes negligible as the mean increases by the strong law of large numbers. This suggests thatZ n (t) converges to a deterministic vector processz(t) defined bȳ We are now ready to state our theorem for case (iv) in Theorem 3.1, the situation where the epidemic initially grows at a higher rate than the population.
Proof The proof follows the same ideas as Theorem 5.2 in Andersson and Britton (2000) originating from Ethier and Kurtz (2005), so we leave out some details. We also leave out the error term stemming from the approximation N(u)/n(u) ≈ 1. Defineβ = supz ∈[0,1] 3 β (z) which is finite. Further, let M < ∞ satisfy | (β (x) − β (y))| ≤ M|x − y| which is true since all β (·) are differentiable. From Eqs. (4) and (5) it follows that From Gronwall's inequality, we then have The exponent is independent of n and both terms within the parenthesis tend to 0 with n which implies that the left-hand side tends to 0, also when taking supremum over finite intervals. This completes the proof.
An alternative way of definingz(t) is in terms of differential equations rather than integral equations. If we differentiate Eq. (5), we get, after a bit of algebra, Writing these term by term results in the following set of equations; This set of differential equations are, in fact, well-known, they correspond to the SIR model with demography in a non-growing population (e.g., Diekmann et al. 2012, Sect. 4.4). The starting point was defined asz 0 = (1 − ε 1 − ε 2 , ε 1 , ε 2 ), where ε 1 and ε 2 were assumed small but otherwise arbitrary. Because our main focus lies in the state of the process after a long time, the endemic level, we look for solutions (ŝ,î,r) where all three derivatives are 0. Equating all derivatives to 0 gives the solutions (1, 0, 0) and (ŝ,î,r), wherê Note that we have assumed that α > λ − μ (super-critical case) ensuring that 0 < s,î,r < 1, and we also have thatŝ +î +r = 1. Now, irrespective of the starting point, as long asī(0) > 0,s(0) ≥ 0,r(0) ≥ 0 ands(0) +ī(0) +r(0) = 1, we have that the disease free solution (1, 0, 0) is unstable, and therefore not a limit point of the process. For relevant parameter values (e.g., that the average infectious period is much shorter than the average life-length) this system (7) is known to exhibit damped oscillations and the process converges to the endemic level: (s(t),ī(t),r(t)) → (ŝ,î,r) as t → ∞, (Diekmann et al. 2012, Sect. 4.4). If we combine our results from the current and previous section for the case α > λ − μ, we have that I n (t) first behaves like I (U ) (t). It may hence die out quickly, which it does with a probability tending to (δ + μ)/κ. With the remaining probability (tending to) 1−(δ +μ)/κ, I n (t) grows at a higher rate than N n (t): I n (t) ∼ e (κ−(δ+μ))t vs N n (t) ∼ ne (λ−μ)t . The sandwich coupling of I n (t) between I (L) (t) and I (U ) (t) breaks down when the number of infected exceeds ε 0 N n (t) ∼ ε 0 ne (λ−μ)t for the first time, which hence equals the first time that I n (t) = ε 0 N n (t) ≥ ε 0 ne (λ−μ)t . From the growth rates just mentioned it can be deduced that this time t n satisfies (Note that the denominator is strictly positive by assumption.) At this time-point we have I n (t n ) ≈ (ε 0 n) (κ−(δ+μ))/(κ−(δ+λ)) .
By our assumptions that λ > μ and κ > δ + λ, it follows that the exponent is larger than 1. So, if ε 1 and ε 2 (the starting point of the approximation of the "proportionprocesses" of infectives and recovered) are chosen small enough, then the second approximation kicks in before the first coupling is broken. As a consequence, we then have convergence to endemicity. We summarize the results in the following theorem: Theorem 3.3 Assume that α > λ − μ, let Z n (t) = (S n (t), I n (t), R n (t)) be the epidemic process starting with (S n (0), I n (0), R n (0)) = (n − 1, 1, 0) and letZ n (t) = Z n (t)/ne (λ−μ)t . Then, for fixed n and as t → ∞ we have that Z n (t) → (∞, 0, 0) (the disease free state in a growing population), or else that Z n (t) → (∞, ∞, ∞). The probability of the first event π n := P (lim t Z n (t) → (∞, 0, 0)) satisfies lim n π n → (δ + μ)/κ.
In the latter event in the theorem, we conjecture thatZ n (t) → (ŝ,î,r) as t → ∞ (the endemic state in a growing population), this is true if the system (7) has no limit cycles in the positive octant.

Remark
The system (7) does not allow for i(t) = 0 if the epidemic takes off. Since n(t) is increasing, this implies that I (t) can with large probability be kept above every desired minimum by choosing n(0) = N 0 (t) large enough. This implies that I n (t) → ∞ as t → ∞ as claimed in part (iv) of Theorem 3.1.

Extension to the Uniform Markovian Dynamic SEIR Epidemic Model
Most infectious diseases have a latent state, i.e., a state where an individual has been infected (Exposed) but is not yet able to spread the disease, before becoming infectious. Accordingly, we now extend our model to an SEIR epidemic allowing for this (the added "E" is for exposed). We have the same population model as before. The only difference is in the epidemic model that an individual who becomes infected is first Exposed (or latent) for an exponentially distributed time with rate ν (unless it dies before this) before turning infectious. After the latent period the individual becomes infectious having infectious contacts just like before. If we let ν → ∞ we hence retrieve the SIR model defined above. The transitions and their rates are given in Table 2.
It is not hard to show that the basic reproduction number, defined as the expected number of infectious contacts an infected individual has, equals The first factor is the probability that the individual does not die during the latent state, and the second is, as before, the expected number of infectious contacts while being infectious. As for the Malthusian parameter α, the exponential rate at which the epidemic initially grows, it is obtained from the defining equation ∞ 0 e −αt c(t) dt = 1, where c(t) is the expected rate at which an infected individual has infectious contacts t time units after he/she was infected. In the SEIR model, this contact rate equals κ, Table 2 The uniform Markovian dynamic SEIR epidemic: type of events, their state change (the old state z = (s, e, i, r) is hence changed to z − ) and their rates for the same n(t) = ne (λ−μ)t as before. The corresponding set of differential equations ares (t) = λ 1 −s(t) − κī(t)s(t), e (t) = κī(t)s(t) − (λ + ν)ē(t), i (t) = νē(t) − (λ + δ)ī(t), r (t) = δī(t) − λr(t).
As in the SIR case, this set of differential equations corresponds to the SEIR model with demography for a non-growing population (Diekmann et al. 2012, Exercise 4.8). The endemic equilibrium (when α > λ−μ) is obtained by setting all derivatives equal to 0 and finding the unique positive solution. It is given by (ŝ,ê,î,r where b = ν/((λ + ν)(λ + δ)). The theorems for the endemic situation apply also to the SEIR case using identical method of proofs, but with the new α from Eq. (12), and for the new vector processZ(t) and its deterministic limitz(t) with its endemic equilibrium just defined. As for the initial phase the limiting probability of a minor outbreak lim n π n is different from the SIR model and equal to (δ + μ)/κ + μ/(ν + μ) as mentioned above.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.