Epidemic modelling by birth-death processes with spatial scaling

In epidemic modeling, interpretation of compartment quantities, such as s , i , and r in relevant equations, is not always straightforward. Ambiguities regarding whether these quantities represent numbers or fractions of individuals in each compartment rise questions about signiﬁcance of the involved parameters. In this paper, we address these challenges by considering a density-dependent epidemic modelling by a birth-death process approach inspired by Kurtz from 1970s’. In contrast to existing literature, which employs population size scaling under constant population condition, we scale with respect to the area. Namely, under the assumption of spatial homogeneity of the population, we consider the quantities of susceptible, infective and recovered per unit area. This spatial scaling allows diﬀusion approximation for birth-death type epidemic models with varying population size. By adopting this approach, we anticipate to contribute to a clear and transparent description of compartment quantities and parameters in epidemic modeling.


Introduction
A classical approach to epidemiological modelling is by the use of ordinary differential equations in which the population is divided into different compartments; see for instance [2] for a historical expose.Those equations attempt to capture empirical dynamics.The scaling is however not always clear; for instance it can be a question of whether the involved variables describe fractions or densities of individuals in the different compartments.Scaling is an inherent problem for ordinary differential equations describing population dynamics or epidemiological evolution of small population sizes whence the individuals are discrete quantities in contrast to the solutions of the ordinary differential equations.
Another inherent problem for ordinary differential description, at least for small populations, is that prediction error is not taken into consideration, whence models that involve randomness is looked for.An often used stochastic approach is to the right-hand side of an epidemiological ordinary equation add a noise representing future uncertainty or a random environment.Such equation is then formulated as a stochastic differential equation.The added noise is then usually of such a type that global existence and uniqueness of a solution can be obtained, see for instance [1,3,4,6].
In [5,[11][12][13][14], the modeling approach is instead based on birth-death type continuous time discrete-valued Markov chains of a certain parameter dependent structure directly applicable to compartment models with a constant total population size where the parameter is the population size.The scaled Markov chains then describe the fraction of individuals in each compartment.For huge population sizes the scaled Markov chains can be approximated by classical epidemiological ordinary differential equations and for intermediate population sizes the scaled Markov chains can be approximated by stochastic differential equations with square root diffusion terms for which a solution of such an equation is defined until the first hitting time of the axes whereafter in [11][12][13][14], the solution of the stochastic differential equation is killed meaning that the stochastic differential equation approximation is not applied after such a hitting time since the epidemiological dynamics may afterwards need another model.By this continuous time Markov-chain approach with constant total population size, the interpretation of the variables of the limiting classical epidemiological ordinary differential equation is clarified as the fraction of individuals in each compartment.
In this paper the population size, due to deaths or migration, is allowed to change by time which means that the population size can not be used a scaling parameter directly applying results in [5,[11][12][13][14].Instead, it is assumed that the number of individuals in each compartment in a region with a given area is described as birth-death type continuous time discrete-valued Markov chain with the area size serving as a scaling parameter.The scaled Markov chain then describes the areal density of individuals in each compartment.Here the scaled time-continuous Markov chains can be approximated by classical epidemiological ordinary differential equations for large region areas and by approximative stochastic differential equations for intermediate region sizes.It means that here, the involved variables of the limiting classical epidemiological ordinary differential equation describe densities of individuals in each compartment.The interpretation of the involved parameters then differs somewhat from the interpretation of the parameters for a population size scaled time-continuous Markov chain.
In Sect.2, birth-death type continuous-time Markov chain modeling is presented.In Sect.3, a critical part for existence and uniqueness of a solution of the continuoustime Markov chain is pointed out.In Sect.4, an initial reproduction number under the population size scaling with constant population size is compared with an initial reproduction number under area scaling.In Sect.4, convergence of scaled time-continuous Markov-chains to corresponding solutions of ordinary differential equations is motivated.In Sect.6, diffusion approximation for the case with no emigration (but allowed immigration and births and deaths) of area scaled continuous-time Markov chains is considered together with a limiting result of the probability for the diffusion approximation to be outside a box ( , M) 3 , 0 < < M, before a horizon time T. Also in Sect.6, the diffusion approximation is compared to stochastic perturbation of parameters; it means a comparison between an approximative stochastic differential equation of a continuous time Markovchain and a stochastic differential equation for which the noise is added to right-hand side of an ordinary differential equation.In Sect.7, numerical simulations of area-scaled continuous time Markov chains, their diffusion approximations and numerical solution of corresponding ordinary differential equations are illustrated under different values of the parameters describing emigration, immigration, births and deaths.In an additional appendix the situation where all mentioned parameters are included is presented.

Time-continuous Markov-chain epidemiological modelling
A typical introductory example in epidemic modelling is the so-called SIR model, suitable for huge populations, where s(t), i(t), and r(t) are fractions of susceptible, infected and recovered individuals, and the parameters β and γ represent a transmission rate and a recovery rate respectively (see e.g, [2]).Here, the recovered are, for simplicity, supposed to never be infected again.
Applying a general setup in [7,[10][11][12][13][14], a natural model for small populations is, where n is the total population size, s n (t), i n (t), and r n (t) representing again the fractions of susceptible, infected, and recovered.The processes {N -1,1,0 (t)} t≥0 and {N 0,-1,1 (t)} t≥0 are two independent standard Poisson processes, each responsible for one of the two types of discrete events in continuous time, namely, transmissions and recoveries.The term represents the modelled number of infected individuals during the interval [0, t].This can be explained by assuming that the intensity of the number of infected individuals is supposed to be proportional to the population size n and the fractions of susceptible and infected individuals, respectively.One possibility is to identify β as the product of an intensity k of pairwise interactions that a single individual has and a probability p that a susceptible-infective encounter leads to transmission i.e. β = kp.The average number of contacts that an individual has per time unit is then k and the probability that a meeting between a susceptible and an infective results an infection is p.The intensity with which a single susceptible is infected is then k times the empirical probability i n that the other one is infected times the probability p that the other one will infect.This means that the intensity with which a single susceptible is infected is ki n p = βi n .The total intensity of interactions leading to infection is then βi n times the number of susceptible ns n , namely, βi n ns n .The term represents the modelled number of recovered individuals during the interval [0, t].This can be explained by the fact that the intensity at which recoveries occur is proportional to the number ni n of infected individuals.The coefficient γ is then the intensity at which one single individual recovers, i.e., the time to recovery is exponentially distributed with a mean of 1/γ .The system (1) is conservative in the sense that if s(0) + i(0) + r(0) = 1 then indeed s(t) + i(t) + r(t) = 1 for t > 0 which means that for s(t), i(t), and r(t) to be fractions indeed make sense.Similarly the system (2) is conservative in the sense that if s n (0) + i n (0) + r n (0) = 1, then also s n (t) + i n (t) + r n (t) = 1 for t > 0.
For non-conservative SIR models for which s(0) + i(0) + r(0) = 1 may not imply s(t) + i(t) + r(t) = 1 for t > 0, for instance if -= α 1 = α 2 = α 3 in (4) is not valid, then an interpretation of s(t), i(t), and r(t) as fractions are in doubt.If < 0 and α = min(α 1 , α 2 , α 3 ) < 0 then the sum of s(t), i(t), and r(t) is even negative for t > -ln( /( + α))/α.However, if ≥ 0, and s(0), i(0), and r(0) are all positive, by standard arguments s(t), i(t), and r(t) are all positive (for t < τ = inf{u > 0 : Similarly, a finite population discrete event continuous time system such as (2) but with varying population size due to the migration needs some contemplation.Instead of scaling with respect to a varying population size we suggest considering a region with area A and let s A , i A and r A be the densities, i.e., the numbers in the three compartments divided by A.Here it is natural to express the migration rate as the difference rates between immigration and emigration rates = + --and similarly it is natural to express the differences between birth rates and death rates explicitly, where independent Poisson processes N +,mig 1,0,0 (mig as in migration), N -,mig 1,0,0 , N + 1,0,0 , . . .N + 0,0,1 are added, and It means that s A (t), i A (t), and r A (t) are not negative for t < τ A .The term is the number of interactions in the region that leads to infection in the region.The rate with which infection occurs is proportional to the area A, the concentration s A of susceptible, and the concentration i A of infected.This is similar to chemical reaction modelling, but here with areas instead of volumes.The quicker the movements of individuals the bigger β, and the quicker the infection is spread.
A more specific interpretation of β can be to think that around each infective individual there is a region, a 'contagion region' , such that a meeting with a susceptible individual within that region leads to infection with intensity α, i.e. a meeting with a susceptible individual in such a region leads to transmission with probability αdt within an infinitesmal time period of length dt.With a the mean area of such a contagion region we let β = aα.With s A the number of susceptible per area, within the 'contagion region' the expected number of interactions with a susceptible individual is as A .The intensity of infections that an infective infects is as A α = βs A .With Ai A the total number of infective, the total intensity of infections is βs A Ai A .
Parameter β can then be seen as a time-space intensity of interactions that lead to infection.
Each of the other terms can be explained separately, for instance, is the number of immigrated individuals to the region in [0, t]; + is the immigration rate per area unit and A + is the immigration rate to the region.Also for instance the term is the number of born susceptible in [0, t] with each susceptible giving births with intensity λ 1 , As A (u) the number of susceptible in the region at time u and Aλ 1 s A (u) the total birth intensity of those individuals at time u.The term N - 1,0,0 t 0 Aμ 1 s A (u)du is likewise the number of dead susceptible in [0, t] with each susceptible having death intensity μ 1 , As A (u) the number of susceptible in the region at time u and Aμ 1 s A (u) the total death intensity of those individuals at time u.
For t > τ A , the dynamics have to be reformulated since we should obviously have min(s A (t), i A (t), r A (t)) ≥ 0. The critical term is, for -> 0, N -,mig 1,0,0 t 0 A -du for which the intensity is positive even though s A (t) may be zero.The other terms will have zero intensity at time t if s A (t), s A (t), or s A (t) are zero.

Existence and uniqueness of a solution to (5)
Existence and uniqueness of a solution to ( 5) is shown iteratively, first up to the first jump time given the involved Poisson processes.Until the first jump time the process is constant.Then we get the solution up to the second jump time, observing again that the process is constant between the jump times.It is critical to show that it will not be infinite intensity in finite time.That can be secured by Gronwall's inequality.By adding the lines in model ( 5), with n A (t) = s A (t) + i A (t) + i A (t) where all the terms are non-negative for t < τ A , we get which by a stochastic version of Gronwall's inequality, see e.g., [13,Appendix] implies that sup t<τ A E[n A (t)]e -αt is finite for some α > 0 which in its turn implies finite total intensity of jumps for n A (t) and therefore for s A (t), i A (t), and r A (t).For more general existence and uniqueness result of equations of the type (5), see for instance [7, Chap.6, Thorem 4.1].

R 0 under the population size scaling and density scaling
We consider the number K of individuals that a single initial infected infects where all others are initially susceptible under the population size scaling (2) and the density scaling (5), respectively.We condition on the time τ until recovery of the first individual being exponential distributed with mean 1/γ .

Population size scaling
For the model (2), the intensity with which that single individual infects is the number of pairwise interactions times the probability that the other individual is susceptible times the empirical probability that the other met individual is susceptible, i.e. kps n = βs n .Since initially only one individual is infected, for n large as for the standard R 0 for (1) obtained by the next generation method [15].Knowing that the average time until recovery is 1/γ it means that the mean intensity of infections that the first individual creates is β = kp, the product of the mean number of interactions the first individual has and the probability that such a meeting leads to a new infective.

Spatial scaling
For the model ( 5), the intensity of infections that a single infective infects where all others are susceptible is βs A where we here recall that β = aα, where a is the mean 'contagion region' for an infective and α is the intensity for which a meeting that an infective has with a susceptible resulting into a new infected.For a comparison with (2) we consider the case + = -= λ i = μ i = 0 for i = 1, 2, 3. Then during the time τ , the number of infections that the first infective infects is, also here approximately Po( τ 0 βs A (t)dt).Here As A (0) is the initial number of susceptible in the region with area A. Then During the infection time with mean 1/γ the intensity with which that single infective infects is hence approximately s A (0)β = as A (0)α, where as A (0) is approximately the approximate mean number of susceptible in a contagion region around the infective of which those susceptible will be infected with intensity α.If the initial susceptible are initially uniformly spread with intensity of ρ individuals per area unit, i.e. the number of the susceptible in a region with area A is Po(ρA), then the intensity with which that single infective infects is approximately ρβ = aρα, the mean contagion region times the density times the intensity that a meeting with a susceptible leads to infection.During the infection time with mean 1/τ the mean number of infected is then E[K] ≈ ρβ/γ .

ODE approximation
Similar to (3), we can, by applying the same references as for the derivation of (3), under the assumption for the case -= 0, derive where (s A , i A , r A ) and (s, i, r) are the solutions of ( 5) and (4), respectively.A limiting result similar to (7) for the case -> 0 needs to be modified including stopping times of the type τ = inf{t > 0 : s(t) < } for > 0 as well as a corresponding stopping time for (5).A limiting result for the case -> 0 is therefore beyond the scope of this paper.
To get an intuitive understanding of (7), we can make a compensated Poisson representation of (5).In order to get not too long formulae, we consider the case for which also + = -= λ i = μ i = 0 for i = 1, 2, 3.The proof is almost identical for the general case, and can be found in the Appendix.With the compensated standard Poisson processes we can then write (1) as By Doob's martingale inequality, for each finite T > 0, P-a.s.

Diffusion approximation
We only consider diffusion approximation of ( 5) under the assumption -= 0 since with -> 0 an analysis with stopping times as for the ODE approximation in the previous section is needed which is beyond the scope of this paper.Also here we study for convenience the case for which also + = λ i = μ i = 0 for i = 1, 2, 3.The case for which some of those parameters + , λ i , μ i are non-zero is also here almost identical.We scale the compensated Poisson process with respect to the area appropriately, Donsker's functional central limit theorem says that for large A, and where W -1,1,0 and W 0,-1,1 are independent standard one-dimensional Brownian motions, and the approximation is uniform on bounded time intervals.This inspires to a diffusion approximation {x Note that existence, uniqueness and convergence as A → ∞ in (10) are not trivial.Particularly, the existence of a solution is not evident due to the non-standard nature of equation (10).
For existence and uniqueness of a solution of (10), we first consider the following SDE (which will later be translated into (10)): where B -1,1,0 and B 0,-1,1 are independent standard one-dimensional Brownian motions.
For given M > 0, existence of a solution of ( 11) is obtained until the first exit time of (0, M) 3 under which the drift and dispersion are bounded and continuous, see for example [8, Theorem 2.2, p. 155].For given ∈ (0, M) uniqueness of a solution of ( 11) is obtained until the first exit time of ( , M) 3 for which the drift and dispersion are Lipschitz continuous, see for example [8, Theorem 2.2, p. 155].
For convergence as A → ∞, a procedure can be as follows.First a probability space with two Brownian motions B -1,1,0 (t) and B 0,-1,1 (t) are given.Then the Brownian motions W -1,1,0 (t) and W 0,-1,1 (t) are constructed as above.Then it is possible to extend the probability space so that two independent standard Poisson processes N -1,1,0 (t) and N 0,-1,1 (t) can be defined such that for the compensators Ñ-1,1,0 (t) = N -1,1,0 (t)t and Ñ0,-1,1 (t) = N 0,-1,1 (t)t, have finite exponential moment (see e.g.[7, Corollary 5.5, p. 359] and [12,Lemma 3.12] with references emanating from [9]; see also [14, p. 237] for a somewhat different formulation on the difference between constructed compensated Poisson processes and Brownian motions), resulting into, for log(A) > 1 and T > 0 and x A (0) = xA (0), there exists a random variable K such that where τ A is the first exit time of x from ( , M) 3 and E[exp(λK)] < ∞ for some λ > 0, see e.g.[12, Lemma 3.7 and Theorem 3.13] and [14, Theorem 1.1 and Theorem 1.2], where in the last reference, the box ( , M) 3 is replaced by a region which is allowed to depend on the scaling parameter which here is the A.
The exit probability of xA from ( , M) 3 in [0, T] is small if A is large.We have for instance: For a fixed horizon T , where In particular, ]dt + dM.
We get We also have Hence The diffusion approximation of (5) allowing also -to be positive is where B -1,1,0 , . . .B - 0,0,1 are independent standard one-dimensional Brownian motions.
Remark 1 (Generator argument of ODE-and SDE approximations) Heuristic generator arguments why ODE-approximation and SDE-approximation can be used for the Poissonian SIR model ( 5) can be as follows.We consider for simplicity the case for which + = -= γ i = μ i for i = 1, . . ., 3. The critical case -> 0, not taken care of here, needs a finer reasoning.In terms of transition probabilities we have for the case + = -= γ i = μ i for i = 1, . . ., 3: For f ∈ C 2 (R 3 ; R), ∇f its gradient and H its Hessian we then get for If we neglect terms of order 1/A or higher then at the grid points, it is the generator of ( 1).
If we neglect the o(1/A) remainder term, then, at the grid points, it is the generator of ( 11) up to the first hitting time τ of the axes: τ = inf{t > 0 : s(t)i(t)r(t) = 0}.
Remark 2 (Stochastically perturbed SIR modeling: parameter perturbation) There is a rich flora of stochastic models that are based on deterministic ODE models for which some parameters are subjected to noise, we mention among others, [1,3,4,6].
For instance consider the SIR ODE (1).If say β is not precisely not known and for simplicity is randomly fluctuating around a constant β 0 , for example β is replaced by β 0 + σ dW dt where dW dt can be interpreted as a Gaussian white noise i.e. βdt is replaced by β 0 dt + σ dW , where W is a Brownian motion, we then get the SDE In (17) there are no square root dispersion terms as opposed to (11).Starting in (0, ∞) 3 , the solution of (17) remains in (0, ∞) 3 , see for instance [1,3,4,6], while for (11), the boundary of (0, ∞) 3 may be hit.
For the SDE (17), the uncertainty is added to the ODE (1).For the SDE (11) the uncertainty is instead built into the model taking into account that spread of diseases is in its nature uncertain.

Simulations
Simulations of the continuous-time Markov chain (5), the SDE (16), and the ODE (4) for some different values of + , -, λ i , and μ i for i = 1, 2, 3, are here presented.The simulations of (5) are exact in the sense that there are no numerical errors.For the SDE (16) and the ODE (4), Euler (Maruyama) approximations are applied.
For the case -> 0, the simulations are stopped just before the susceptible part is negative wherafter ( 5), (16), and (4) are not valid anymore.
For the case -= 0, non-negativity for ( 5) and ( 4) do not occur while for (16) a solution is defined until the first hitting time of the axes.The behaviour of the solution of ( 16) near the boundary of [0, ∞) 3 is not evident.In the simulations it was found reasonable to let the numerical solution continue after hitting the axes by applying the x → max(x, 0) function at each Euler Maruyama step, i.e. a projected version of the Euler Maruyama scheme is applied.It means that the numerical solution solves a constrained version of the SDE (16) inwards reflecting at the boundary of [0, ∞) 3 .It is then not obvious whether the reflected version of (16) need additional 'push-terms' (described by the local times at the boundary) or is instantaneously reflecting without push-terms.An alternative simulation could, for the case when also + = 0, to let {0} be absorbing for both the susceptible and the infective compartments while allowing {0} to be reflecting for the recovered compartment allowing the dynamics to evolve even if initial value for the recovery compartment is zero.A third alternative could be to allow the numerical simulation to take negative values and replacing the square roots by square roots of absolute values.The SDE ( 16) is after all not a model but an approximation of (5).
In all the simulations (see Fig.

Conclusion
To conclude, in this study we acknowledge that the population size changes over time due to factors like births, deaths or migration.This means that using population size as a scaling parameter, as follows directly from previous studies [11][12][13][14], is not suitable.Instead, we propose a different approach: we model the number of individuals in each compartment in a region of a specific area using a birth-death type continuous time Markov chain.In this setup, the area size becomes the scaling parameter, and the resulting scaled Markov chain shows the density of individuals in each compartment.This approach allows us to approximate the scaled time-continuous Markov chains using classical ordinary differential equations for large regions, and approximating stochastic differential equations for intermediate size regions.We underline that for the case with constant total population size, the variables in a classical epidemiological ordinary differ-  In summary, our research provides a detailed analysis of epidemic modeling techniques, focusing on understanding compartmental quantities like susceptible s, infected i, and recovered r.We introduce a new density-dependent approach inspired by Kurtz's from 1970s' to clarify the significance of key parameters.Unlike traditional methods that rely on constant population size scaling, our approach uses area-based scaling, enabling us to adapt to fluctuations in population size.Through this method, we aim to provide a clearer and more transparent illustration of compartment quantities and parameters in epidemic modeling.
Potential future directions involve investigation of the distributions of the first hitting times for the approximate diffusion process on the axes as well as investigation of continuations of solutions after the first hitting times and, for the case with emigration, the first hitting time of the axes for the continuous-time Markov chain.

Figure 6
Figure 6 50 Euler-Maruyama approximative simulations of spatial scaling of the SDE approximation (16) together with the ODE approximation (4) respectively, indicating the probability density at each time point for three cases: (a) + = -= λ i = μ i = 0 for i = 1, 2, 3; (b) μ 1 = 4, + = -= λ 1 = μ 2 = λ 2 = μ 3 = λ 3 = 0; (c) -= 2, + = λ i = μ i = 0 for i = 1, 2, 3.The scaling of the axes in the figures are also the same.The figures here are similar to those in Fig. 5.In (c), the simulations are stopped just before the susceptible part is negative.The probability density of that exit time is also here indicated.The infective component and the recovery component are not stopped for the Euler-Maruyama approximations; instead they are conveniently projected in order for them to be non-zero similar to the continuous-time Markov chain; the dynamics stops first at the exit time for the density of susceptibles