A Novel Stochastic Epidemic Model with Application to COVID-19

In this paper we propose a novel SEIR stochastic epidemic model. A distinguishing feature of this new model is that it allows us to consider a set up under general latency and infectious period distributions. To some extent, queuing systems with infinitely many servers and a Markov chain with time-varying transition rate are the very technical underpinning of the paper. Although more general, the Markov chain is as tractable as previous models for exponentially distributed latency and infection periods. It is also significantly simpler and more tractable than semi-Markov models with a similar level of generality. Based on the notion of stochastic stability, we derive a sufficient condition for a shrinking epidemic in terms of the queuing system's occupation rate that drives the dynamics. Relying on this condition, we propose a class of ad-hoc stabilising mitigation strategies that seek to keep a balanced occupation rate after a prescribed mitigation-free period. We validate the approach in the light of recent data on the COVID-19 epidemic and assess the effect of different stabilising strategies. The results suggest that it is possible to curb the epidemic with various occupation rate levels, as long as the mitigation is not excessively procrastinated.

more general, the Markov chain is as tractable as previous models for exponentially distributed latency and infection periods.It is also significantly simpler and more tractable than semi-Markov models with a similar level of generality.Based on the notion of stochastic stability, we derive a sufficient condition for a shrinking epidemic in terms of the queuing system's occupation rate that drives the dynamics.Relying on this condition, we propose a class of ad-hoc stabilising mitigation strategies that seek to keep a balanced occupation rate 1. Introduction Classical epidemic models (e.g., Ross, 1916;Kermack et al., 1927;Hethcote, 2000) are powerful tools to understand the spread of diseases and support public health policies.However, these models are deterministic and consequently do not capture the underlying uncertainties of the spread.Stochastic epidemic models (Allen, 2008;Britton, 2010) were therefore designed to better capture some of these uncertainties and provide more realistic support for decision making.
Perhaps due to its simplicity, the susceptible, infected, removed (SIR) is arguably the most utilised class of stochastic epidemic models.Trapman and Bootsma (2009) made use of this framework to demonstrate the advantage of an M/G/1 queuing model to estimate the size of an epidemic at the time of detection.Some time later, classical M/M/S queues were utilised to estimate the whole outbreak of the Ebola virus (Dike et al., 2016).More recently, Barraza et al. (2020) employed pure Birth processes to fit data from the initial stages of the COVID-19 epidemic.Underpinning the analysis is the theory of Markov processes (Brémaud, 1999), used to model the transition of individuals among the SIR populations.
Markov models provide a powerful analytical framework for SIR models, allowing for example the treatment of non-homogeneous populations (López-García, 2016).One of the limitations, however, is that the duration of the infectious period is assumed exponential, thus narrowing the technique's applicability (Clancy, 2014;Gómez-Corral and López-García, 2017).A possible alternative is to develop more complex block-structured Markov chains that can mimic certain types of non-exponential infection times (Lefèvre and Simon, 2020;İşlier et al., 2020).Albeit limited, the added flexibility comes at the price of less interpretable and tractable models.
A thorough treatment of general infection times demands a class of semi-Markov processes known as piecewise-deterministic processes (PDP) (Davis, 1993).Clancy (2014) uses such a class and martingale theory to derive the distribution of the number of infected individuals throughout the pandemic under generally distributed infection times.Later, Gómez-Corral and López-García (2017) utilised a similar model to obtain the distribution of the number of secondary cases due to an infected individual.However, it requires a memory of the disease progression of all currently infected individuals, which impacts the model's analytical and computational tractability and limits its use.
Inspired by the recent COVID-19 outbreak, this paper uses the SEIR (susceptible, exposed, infected, removed ) framework, which is analogous to SIR except for considering a latency period, whereby a recently infected individual does not manifest the disease nor can transmit it for a limited period after acquiring the illness.In addition to being more general, such a modelling choice is adequate for the COVID-19 epidemic, given its non-negligible latency period (Backer et al., 2020).Examples of stochastic SEIR models (Artalejo et al., 2015;Lopez-Herrero, 2017;Amador and Lopez-Herrero, 2018) assume exponentially distributed latency and infection periods.The first work concerns the duration of the outbreak, the second the transmission rate per infected person and the latter the epidemic size.This paper builds upon previous literature (Trapman and Bootsma, 2009) with innovative use of two M t /G/∞ models to describe the epidemic's stochastic behaviour.However, whilst the M/G/1 queue in (Trapman and Bootsma, 2009) models the epidemic up to the time of detection, the proposed model covers the whole outbreak.Such a generalisation hinges on two unique novelties.Firstly, we describe the input process as a time-varying Poisson process, which enables us to describe the variation of infection rates as a function of the system's dynamics.Secondly, by realising that there are no limits in the number of new infections and considering that conditions progress in parallel, we select a queuing model with infinitely many servers.
To the best of our knowledge, this is the first work to introduce a timevarying Markov chain capable of emulating a stochastic epidemic model with general latency and infection times.This innovation relies on the results of Eick et al. (1993), that demonstrate that the output process of M t /G/∞ queues comprises a series of time indexed Poisson variables.Albeit as tractable as the models that assume exponential (memoryless) infection and latency periods (e.g., Artalejo et al., 2015;Lopez-Herrero, 2017;Amador and Lopez-Herrero, 2018), the proposed framework is more general than the complex block-structured models (e.g., Lefèvre and Simon, 2020;İşlier et al., 2020) as it does not impose any assumption on the infection and latency period distributions.
Even if our approach's level of generality is matched in part by previous PDP models within the more restricted SIR framework (Clancy, 2014;Gómez-Corral and López-García, 2017), these require memory of the disease progression of all individuals in the infected population.Of course, this is infeasible for all but tiny population sizes and limits such models' applicability to support decision making.In contrast, the proposed framework considers general latency and infection periods by keeping track solely of the new expositions within a complete cycle of the disease: from catching the virus to entering the removed population.
In addition to the methodological innovations, we propose a novel strategy to curb the epidemic that is based on the classical occupation rate of the M t /G/∞ exposition-to-removal queue.We claim that this measure is analogous to the deterministic reproduction number as it also provides a dynamic estimate of the epidemic's short-term growth.The strategy relies on mitigating actions specially tailored to maintain the occupation rate ρ of the system as close as possible to a prescribed level 0 < ρ < 1 at all times and thus ensure a shrinking epidemic.
We use COVID-19 data from the literature to validate the strategy by means of experiments designed to illustrate the utility of the approach whilst providing invaluable insights into the system's response to the proposed class of mitigation policies.The results demonstrate that by maintaining adequate occupation level targets, the epidemic can be conquered, as long as the mitigation is not excessively delayed.This work is organised as follows.Section 2 introduces the mathematical formulation of the proposed SEIR model, starting with the proposed M t /G/∞ queues.We then use these queues' resulting input and output processes to propose a continuous time Markov chain to describe the epidemic's evolution.
Section 3 analyses stochastic stability and makes use of the results to propose a class of ad-hoc mitigation strategies to curb the epidemic.Section 4 employs COVID-19 data from the literature to establish a set of experiments designed to illustrate the performance of the model in a real-world setting.The experiments also highlight the effectiveness of prescribed mitigation policies belonging to the class introduced in the previous Section.Finally, Section 5 concludes the manuscript.

A Few Notations
Throughout this paper we shall be using the following notation.Let Z + be the usual set of non-negative integers and N = {0, 1, 2, 3, ...N } ⊂ Z + , where N is a finite number.Consider Ω = N 4 and define the probability space (Ω, F , P).
In addition, E(•) stands for the mathematical expectation.

Mathematical formulation
This section proposes a stochastic dynamic formulation in a probability space (Ω, F , P).It relies on the classical deterministic SEIR epidemiological model (Ross, 1916;Kermack et al., 1927;Hethcote, 2000), thus depicted in Figure 1 and in Table 1.
The model comprises four compartments, namely: (i) susceptible, (ii) exposed, (iii) infected and (iv) removed.Susceptible individuals can be infected if In the deterministic model in Figure 1, susceptible individuals become ill at a rate β > 0 upon making contact with infected (infectious) individuals, thus resulting in an infection rate δ = βSI.The individuals who acquire the condition immediately become exposed and initiate their latency period, which lasts on average σ units of time.Upon completing the latency period, exposed individuals become infectious for an average of γ time units.After that, they enter the removed population either due to acquired immunity or death.Table 1 below presents the parameters of the deterministic model:

Stochastic formulation
Although the model in Figure 1 is invaluable to understand the underlying process, the dynamics are truly stochastic.Both the latency and the recovery period are, indeed, stochastic variables.See for example Backer et al. (2020) and Verity et al. (2020) for estimations of the latency and recovery periods at the early stages of the COVID-19 epidemic.In this Section, use the insights of Trapman and Bootsma (2009) to model the SEIR dynamics as a queue with infinite service capacity.
Queue 1 in Figure 2 represents an individual's trajectory from acquiring the disease and therefore becoming exposed, to their removal of the system.
Because the individuals' trajectories are assumed independent, the system can be modelled as a M t /G/∞ queue (Eick et al., 1993), since there is no upper limit on the number of simultaneous infections.The service time is the sum of two generally distributed random variables σ and γ, representing the latency and recovery periods, see Table 1; the latter corresponds to the length of the infectious period.

The arrival and departure processes
We start by modelling the arrival process shared by both queues.At any given time t ≥ 0, the rate of contagion (exposition) is given by where β is a scalar parameter, S(t) and I(t) are the sizes of the susceptible and infected populations at time t, respectively.
By design, the output of Queue 2 -in Figure 2 -is the number of new infections, i.e. individuals become infected when they depart this queue.Let δ e (t) be a random variable denoting the number of departures from Queue 2 at time t ≥ 0. According to Eick et al. (1993), δ e (t) is a Poisson random variable with mean: This is because the latency period σ is also the service time of Queue 2.
Let δ i (t) be a random variable representing the number of departures from Queue 1 at time t ≥ 0. Each departure is a newly removed individual and the total time from exposition to removal is the sum of random variables σ and γ, which is also the service time of Queue 1. Once again, the results of Eick et al. (1993) imply that δ i (t) is a Poisson random variable with mean: (3)

Markov formulation with time-varying transition rate
We make use of the arrival and departure rates of Queues 1 and 2 in Eq.
(1)-( 3) to define a time-varying Markov process X t , t ≥ 0 that describes the evolution of the populations in the SEIR compartments.At any time t ≥ 0, By definition, X t , t ≥ 0 will be subject to random jumps that happen whenever either of these events occur: i) a new exposition, ii) a new departure from Queue 2 or iii) a new departure from Queue 1.
To describe the dynamics of the process X t , t ≥ 0, we make use of the fact that the input and output of both queues are described by Poisson variables (Eick et al., 1993) at any given time t ≥ 0. Consequently, the time until the next event will be exponentially distributed and the total jump rate at time t ≥ 0 is: Let {τ 0 , τ 1 , . ..} be the sequence of jumps in the system, with τ 0 ≡ 0 and τ k+1 > τ k , ∀k ≥ 0. Now, assume a jump occurs at time t = τ k .Then, the following holds: (5) The first expression on the right-hand side of Eq. ( 5) corresponds to a new contagion/exposition, which happens at rate λ(t), see Eq. ( 1), and implies the transference of an individual from the susceptible to the exposed population.
The second expression corresponds to a new infection that happens upon the departure of an individual from Queue 2, which occurs at rate δ e (t), see Eq. ( 2).
In that case, this individual moves from the exposed to the infected population.
Finally, the third possibility is the departure of an individual from Queue 1, which happens at rate δ i (t), see Eq. ( 3).In that case, this individual migrates from the infected to the removed population.
Finally, after the jump at t = τ k , k ≥ 0, the value of the exposition rate λ(t) also changes, and becomes: where the new values on the right-hand side of the expression above vary as a function of the transition probabilities in (5).Clearly, λ(t) remains unaltered between successive jumps.Given an initial distribution and the transition probability in (5), the process X t , t ≥ 0 is a continuous-time Markov chain.

Stochastic stability and reproduction number
Whereas deterministic models rely on the evaluation of trivial equilibrium to derive stability conditions and the so-called reproduction number R 0 (e.g. van den Driessche and Watmough, 2002), the notion of stochastic stability (e.g., Meyn and Tweedie, 1993) provides an ideal framework to evaluate the conditions for a receding epidemic as time elapses.Queuing theory connects this notion with the so-called occupation rate, a measure of the input to output ratio as time elapses (Shortle et al., 2018).
Consider Queue 1 in Figure 2. Assuming a very large population, the stability of such a queue hinges on the occupation rate (e.g., Shortle et al., 2018): and can be ascertained if a finite time t ≥ 0 exists such that ρ(t) < 1, ∀t ≥ t.
This guarantees that the system stabilises and the number of customers in the queue remains finite.With a finite population, however, stability is guaranteed because the number of infected individuals will remain within finite, albeit possibly large, bounds.In that case, we are interested in the trend of Queue 1, which indicates whether the epidemic is increasing or receding.Theorem 1 below establishes the condition for a receding epidemic.
Theorem 1.Consider the Markov process described in Section 2.1.2and assume that ρ(t) < 1 for all t > t ≥ 0, with t < ∞.Then, it follows that: where {τ 0 , τ 1 , . ..} is the sequence of jumps in the system, as defined in Section 2.1.2.
As a consequence of Theorem 1, we can interpret ρ(t) as the reproduction number of the system at time t and ρ(t) < 1 ∀t > t ≥ 0 as a sufficient condition for the epidemic to stabilise and shrink.It is worth pointing out that the theorem indicates that the epidemic is receding only when the sum of exposed and infected individuals, i.e. all latent and manifested infections, decreases.It highlights that, although important, the number of infected individuals is not sufficient to establish the tendency of the epidemic.Indeed, it is quite intuitive that the latter depends upon the number of latent infections, which are bound to occur shortly and therefore cannot be ignored.
The model indicates that we have to consider the disease's whole cycle, from exposition to removal, to evaluate the epidemic trend, such as in Eq. ( 7).
It underscores the importance of keeping an accurate track of the epidemic's evolution through an efficient testing strategy.The accurate evaluation of such a cycle's length will also be essential to evaluate the lags between mitigating actions and decreases in the number of infections and expositions.Although mitigation can prevent future expositions and infections, it has no effect on recent transmissions still in the latent stage, which will continue to manifest and may actually drive infection up in the first stages of the mitigation.This will become more evident in the experiments carried out in Section 4.

Mitigation strategies
Theorem 1 provides a basis for developing mitigation strategies to stabilise the epidemic, whereas the Markov model in Section 2.1.2enables us to evaluate the long-term effects of such strategies.Following the literature, we introduce the control (mitigating action) in the form of non-pharmaceutical interventions to limit the spread of the disease (e.g., Ferguson et al., 2020;Kantner and Koprucki, 2020;Tarrataca et al., 2021).A control level 0 ≤ u(t) ≤ 1 attains a reduction of 100u(t) percentage points in the transmission rate at time t ≥ 0, thus resulting in a controlled exposition rate: Let π = {u(t), t ≥ 0} be a mitigation strategy and let Π denote the set of all feasible strategies Π = {π : 0 ≤ u(t) ≤ 1, ∀t ≥ 0}.The dynamics of the system under any strategy π ∈ Π can be evaluated by running the Markov chain X t , t ≥ 0 with the same underlying dynamics, but making λ(t) = λ(t, u(t)), t ≥ 0.
We are particularly interested in the class of stabilising policies Π S ∈ Π, such that: where ρ(t) is the uncontrolled occupation rate in Eq. ( 7).It is easy to see that such policies lead to a controlled input rate ρ c (t) = λ(t,u(t)) δi(t) ≤ ρ, ∀t > t, see Eq. ( 9).Henceforth, ρ will be called target occupation level and π = ( t, ρ) will describe any stabilising policy π ∈ Π S that satisfies Eq. ( 10).In the next Section we will evaluate these policies in the light of the COVID-19 epidemic and explore the stabilising effect of parameters ρ and t in Eq. ( 10).Observe that, because the population is finite, the effect of the control is highly sensitive to and can be limited by a delayed start of mitigation.
Observe, however, that the Markov model in Section 2.1.2is not limited to the proposed class of mitigation strategies.In fact, one can substitute any arbitrary mitigation policy π ∈ Π in Eq. ( 9) and run the controlled Markov chain to evaluate the effect of such policy.To demonstrate that, our experiments will also include the on-off lock-down policies proposed by Tarrataca et al. (2021).
Designed to control hospital bed's occupation, these policies trigger a full scale lock-down when infections surpass a prescribed upper bound; conversely, mitigating actions are lifted when infections fall bellow a prescribed lower bound.

Numerical Experiments
This Section validates the proposed stochastic SEIR dynamic model in Section 2.1.2 in the light of data from the COVID-19 epidemic.We replicate the parameters of Tarrataca et al. (2021) concerning the spread of the epidemic in Brazil.The parameters are listed in Table 2. To complete the necessary data to run the model, we also need the distributions of the latency period σ and of the latency plus infection period (σ + γ), which represents the full disease cycle.The distributions are compatible with previous observations (Backer et al., 2020;Verity et al., 2020) and assumed discrete to cope with the typical daily collection of data.Tables 3 and 4 unveil the distributions of σ and (σ + γ), respectively.

The effect of delayed mitigation
The first series of experiments examine policies π = ( t, ρ) ∈ Π S that follow Eq. ( 10) with a fixed target occupation level ρ = 0.95.The objective is to evaluate the effect of delaying the start of mitigating actions by varying the first parameter.Cases A-G in Table 5 feature different triggering times for the mitigating actions, starting with the case where no mitigation is enforced ( t = ∞).For each experiment, we ran a full realisation of the system for two years, starting from the outset of the epidemic.The difference with respect to Case B is a rather steep increase in the final number of that catch the disease at some point, represented by the removed population.This number increases from 0.4% in Case B to around 17% in Case C, thus highlighting the sensitivity of the spread with respect to mitigation delay.The results also show that a three-week delay produces a large increment in the peak of infections, for even though this peak only reaches around 1% of the population, it is about four times as large as it would be if the mitigation strategy started 21 days earlier.

The effect of distinct occupation rates
This Section investigates the effect of distinct target occupation rates ρ in the stabilising policies.To provide an interesting baseline, we start our mitigating actions on t = 105, when the epidemic reaches a significant proportion of the population but can still be controlled as in Case E of Figure 7.The set of experiments in Table 6 assesses the effect of varying the target occupation level from 60% to 90% in regular intervals.To account for the stochastic fluctuations and prevent biases, Figures 10 to 13 show the median of 100 realisations of the stochastic system within a two-year interval.Figure 10 presents the results for Case H.As one can observe, the number of infected people peaks around 18% on day 114.As expected, the stronger mitigation has no effect in the early stages and is not able to prevent a similar peak as in the baseline (Case E).However, as time elapses we can notice that the increased mitigation produces a much more pronounced decrease in infections.
In effect, infections drop to about 2% by day 200 and 1% on the 231 th day.The increased control is able to keep the removed population at about 52% at the end of two years, whereas that level was 97.5% in the baseline.It is noteworthy that the mitigation is kept at a high level after the system stabilises to avoid The trend of increasing infection levels gathers speed in Case J, as the decrease in the required mitigating actions leads to about 93% of removals within the two-year horizon.On day 200, the number of infections reaches 5% of the population, and 117 days later it reaches 1%.Coupled with the significant decrease in the susceptible population, which reduces the potential spread, the decrease in ρ produces a steep decrease in the control levels.These reach a minimum of 5% as the susceptible population starts to stabilise.However, increased levels of mitigation are required later to avoid a second wave.
Finally, Figure 13 presents the results for Case K.As the target occupation level approaches that of the baseline (Case E), the overall behaviour becomes quite similar.We can notice a small decrease in the overall removals, that reach about 3.5% as opposed to 2.5% in the baseline.Despite dropping mitigating

The effect of on-off lock-down policies
This Section briefly evaluates the effect of on-off lock-down policies.Proposed by Tarrataca et al. (2021), these policies demand a full scale lock-down (u(t) = 1, ρ = 0) when the number of infections surpass a prescribed upper limit; conversely, all measures are lifted (u(t) = 0, ρ = ∞) as soon as these numbers drop below a lower bound.We simulate the two policies whose upper and lower bounds appear in Table 7.By drawing a parallel between stochastic stability and the traditional reproduction number, the paper introduced stabilising strategies that can curb an epidemic under very general conditions, provided that the mitigation is initiated in a timely manner.Beyond the academic contributions, we argue that the generality of the proposed approach renders it invaluable to support real-world decision making in the face of future epidemics.The methodology was validated using COVID-19 data from the literature.The results provide a panorama of insights into the pros and cons of distinct stabilising mitigation strategies over a two-year horizon.
As expected, the experiments illustrate that early intervention is vital to prevent the disease from affecting a large portion of the population.However, to effectively prevent the spread, we need very high levels of mitigation over the whole two-year horizon.The results also suggest that delayed mitigation leads to a dramatic increase in the overall number of infections.In effect, delays tend to produce persistent infection levels over time, even in the presence of mitigation, thereby increasing the spread even though the infection curves are effectively flattened.Beyond flattening the curves, this behaviour suggests that we also need to ensure, by acting swiftly, that their summit is tolerable from both societal and health care perspectives.
The proposed approach leaves many research avenues to be explored in future works.One obvious route is to pursue stochastic optimal control policies that somehow address the compromise between health care aspects, societal issues and the economic burden of mitigation strategies.The challenges involve finding a meaningful trade-off among the different elements that decision-makers need to consider, as well as proposing effective formulations that avoid the curse of dimensionality (Powell, 2011(Powell, , 2019) ) to ensure that the problem remains tractable.
Another branch goes into developing filtering and analytical approaches to estimate a system's parameters considering that the available information is delayed and biased, since dependent on local testing and reporting policies.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Epidemic evolution for Case A: the unmitigated spread

Figures 6 Figure 6 :
Figures 6 features the results for Case D. Despite an extra delay of one week with respect to Case C, it is still possible to flatten the curve and prevent an uncontrollable increase of the epidemic.Nonetheless, the extra week produces a peak of infections five times as large as in Case C, with 5% of people being infected 107 days after the outset of the epidemic.The total removed population skyrockets from about 17% in Case C to about 90% in Case D, thus reinforcing

Figure 7 Figure 7 :
Figure7shows the results for Case E that provide further evidence of the deleterious effect of mitigation delay.A further delay of one week with respect to Case D results in approximately three times as many infected individuals at the peak, that occurs on day 114 and amounts to around 14% of infected individuals.Due to the larger peak, we are able to relax the mitigation earlier and more rapidly, reaching a full relaxation just after seven months.However, since there are still susceptible individuals in the population, mitigation has to re-start from day 396 onwards to avoid a recrudescence of the epidemic.The

Figure 8 :Figure 9 :
Figure 8: Epidemic evolution for Case F

Figure 10 :
Figure 10: Epidemic evolution for Case H

Figure 12 :Figure 13 :
Figure 11: Epidemic evolution for Case I

Figure 14 :Figure 15 :
Figure 14: Epidemic evolution for Case L

Figure 16 :
Figure 16: Time in Lock-down in Cases L and M

Table 1 :
Parameters of SEIR dynamics.

Table 2 :
Parameters of the Simulation.

Table 6 :
Second set of experiments.

Table 7 :
Third set of experiments.