A Novel Stochastic Epidemic Model with Application to COVID-19
Edilson F. Arruda, Rodrigo e Alvim Alexandre, Marcelo D. Fragoso, João B. R. do val, Sinnu S. Thomas
aa r X i v : . [ q - b i o . Q M ] F e b A Novel Stochastic Epidemic Modelwith Application to COVID-19
Edilson F. Arruda a,b , Rodrigo e A. Alexandre b , Marcelo D.Fragoso c , Jo˜ao B. R. do Val d , Sinnu S. Thomas e a Department of Decision Analytics and Risk, Southampton Business School,University of Southampton, 12 University Rd, Southampton SO17 1BJ, UK b Alberto Luiz Coimbra Institute-Graduate School and Research in Engineering,Federal University of Rio de Janeiro. CP 68507, Rio de Janeiro 21941-972, Brasil c National Laboratory for Scientific Computation, Av. Gett´ulio Vargas 333,Quitandinha, Petr´opolis RJ 25651-075, Brasil d School of Electrical Engineering, University of Campinas, Av. Albert Einstein 400,Cidade Universit´aria, Campinas SP 13083-852, Brasil e School of Computer Science and Engineering, Kerala University of Digital Sciences,Innovation and Technology, Technocity, Mangalapuram Thonnakkal POThiruvananthapuram, Kerala-695317, India
Abstract
In this paper we propose a novel SEIR stochastic epidemic model. A distin-guishing feature of this new model is that it allows us to consider a set up undergeneral latency and infectious period distributions. To some extent, queuingsystems with infinitely many servers and a Markov chain with time-varyingtransition rate are the very technical underpinning of the paper. Althoughmore general, the Markov chain is as tractable as previous models for exponen-tially distributed latency and infection periods. It is also significantly simplerand more tractable than semi-Markov models with a similar level of general-ity. Based on the notion of stochastic stability, we derive a sufficient conditionfor a shrinking epidemic in terms of the queuing system’s occupation rate thatdrives the dynamics. Relying on this condition, we propose a class of ad-hocstabilising mitigation strategies that seek to keep a balanced occupation rate ∗ Corresponding author. Tel.: +44 023 8059 7677
Email addresses: [email protected] (Edilson F. Arruda), [email protected] (Rodrigo e A. Alexandre), [email protected] (Marcelo D.Fragoso), [email protected] (Jo˜ao B. R. do Val), [email protected] (SinnuS. Thomas)
Preprint submitted to Elsevier February 17, 2021 fter a prescribed mitigation-free period. We validate the approach in the lightof recent data on the COVID-19 epidemic and assess the effect of different sta-bilising strategies. The results suggest that it is possible to curb the epidemicwith various occupation rate levels, as long as the mitigation is not excessivelyprocrastinated.
Keywords:
OR in health services, Stochastic epidemic models, Markovprocesses, Queuing Theory, Stabilising control.
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 pub-lic health policies. However, these models are deterministic and consequentlydo not capture the underlying uncertainties of the spread. Stochastic epidemicmodels (Allen, 2008; Britton, 2010) were therefore designed to better capturesome of these uncertainties and provide more realistic support for decision mak-ing.Perhaps due to its simplicity, the susceptible, infected, removed (SIR) is ar-guably 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/1queuing 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 wholeoutbreak 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 theCOVID-19 epidemic. Underpinning the analysis is the theory of Markov pro-cesses (Br´emaud, 1999), used to model the transition of individuals among theSIR populations.Markov models provide a powerful analytical framework for SIR models, al-lowing for example the treatment of non-homogeneous populations (L´opez-Garc´ıa,2016). One of the limitations, however, is that the duration of the infectiousperiod is assumed exponential, thus narrowing the technique’s applicability2Clancy, 2014; G´omez-Corral and L´opez-Garc´ıa, 2017). A possible alternative isto develop more complex block-structured Markov chains that can mimic certaintypes of non-exponential infection times (Lef`evre and Simon, 2020; ˙I¸slier et al.,2020). Albeit limited, the added flexibility comes at the price of less inter-pretable 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 dis-tribution of the number of infected individuals throughout the pandemic undergenerally distributed infection times. Later, G´omez-Corral and L´opez-Garc´ıa(2017) utilised a similar model to obtain the distribution of the number of sec-ondary cases due to an infected individual. However, it requires a memory ofthe disease progression of all currently infected individuals, which impacts themodel’s analytical and computational tractability and limits its use.Inspired by the recent COVID-19 outbreak, this paper uses the SEIR ( sus-ceptible, exposed, infected, removed ) framework, which is analogous to SIR ex-cept for considering a latency period, whereby a recently infected individualdoes not manifest the disease nor can transmit it for a limited period after ac-quiring the illness. In addition to being more general, such a modelling choiceis adequate for the COVID-19 epidemic, given its non-negligible latency pe-riod (Backer et al., 2020). Examples of stochastic SEIR models (Artalejo et al.,2015; Lopez-Herrero, 2017; Amador and Lopez-Herrero, 2018) assume exponen-tially distributed latency and infection periods. The first work concerns theduration of the outbreak, the second the transmission rate per infected personand 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 stochasticbehaviour. However, whilst the M/G/ M t /G/ ∞ queuescomprises a series of time indexed Poisson variables. Albeit as tractable as themodels 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-structuredmodels (e.g., Lef`evre and Simon, 2020; ˙I¸slier et al., 2020) as it does not imposeany assumption on the infection and latency period distributions.Even if our approach’s level of generality is matched in part by previ-ous PDP models within the more restricted SIR framework (Clancy, 2014;G´omez-Corral and L´opez-Garc´ıa, 2017), these require memory of the diseaseprogression of all individuals in the infected population. Of course, this is in-feasible for all but tiny population sizes and limits such models’ applicability tosupport decision making. In contrast, the proposed framework considers gen-eral latency and infection periods by keeping track solely of the new expositionswithin a complete cycle of the disease: from catching the virus to entering theremoved population.In addition to the methodological innovations, we propose a novel strategy tocurb 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 thedeterministic reproduction number as it also provides a dynamic estimate of theepidemic’s short-term growth. The strategy relies on mitigating actions speciallytailored to maintain the occupation rate ρ of the system as close as possible toa prescribed level 0 < ¯ ρ < M t /G/ ∞ queues. We then use these queues’ resulting input and output processes topropose a continuous time Markov chain to describe the epidemic’s evolution.Section 3 analyses stochastic stability and makes use of the results to proposea class of ad-hoc mitigation strategies to curb the epidemic. Section 4 employsCOVID-19 data from the literature to establish a set of experiments designed toillustrate the performance of the model in a real-world setting. The experimentsalso highlight the effectiveness of prescribed mitigation policies belonging tothe class introduced in the previous Section. Finally, Section 5 concludes themanuscript. Throughout this paper we shall be using the following notation. Let Z + bethe usual set of non-negative integers and N = { , , , , ...N } ⊂ Z + , where N is a finite number. Consider Ω = N and define the probability space (Ω , F , P ).In addition, E ( · ) stands for the mathematical expectation.
2. 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 1and in Table 1.The model comprises four compartments, namely: ( i ) susceptible, ( ii ) ex-posed, ( iii ) infected and ( iv ) removed. Susceptible individuals can be infected if5
E I RβI σ − γ − Figure 1: The classical SEIR model they come in contact with infected individuals. The exposed population, mean-while, have been infected but the infection is still latent. In that phase, theyare not able to spread the disease and do not manifest any symptom. Oncethe latency period is over, exposed persons become infected and may presentsymptoms. Finally, infected individuals either die as a result of the disease orbecome immune to the disease. In both cases, they migrate to the removed population, which consists of individuals that can no longer contract nor spreadthe disease.In the deterministic model in Figure 1, susceptible individuals become illat a rate β > δ = βSI . The individuals who acquire thecondition immediately become exposed and initiate their latency period, whichlasts on average σ units of time. Upon completing the latency period, exposedindividuals become infectious for an average of γ time units. After that, theyenter the removed population either due to acquired immunity or death. Table1 below presents the parameters of the deterministic model: Table 1: Parameters of SEIR dynamics.
Parameter Description Unit β Transmission rate transmissions/encounter σ Latency period days γ Recovery period days
Although the model in Figure 1 is invaluable to understand the underlyingprocess, the dynamics are truly stochastic. Both the latency and the recovery6eriod are, indeed, stochastic variables. See for example Backer et al. (2020)and Verity et al. (2020) for estimations of the latency and recovery periods atthe early stages of the COVID-19 epidemic. In this Section, use the insightsof Trapman and Bootsma (2009) to model the SEIR dynamics as a queue withinfinite service capacity.
Queue 1 in Figure 2 represents an individual’s trajectory from acquiringthe disease and therefore becoming exposed, to their removal of the system.Because the individuals’ trajectories are assumed independent, the system canbe modelled as a M t /G/ ∞ queue (Eick et al., 1993), since there is no upperlimit on the number of simultaneous infections. The service time is the sum oftwo generally distributed random variables σ and γ , representing the latencyand recovery periods, see Table 1; the latter corresponds to the length of theinfectious period. λ ( t ) δ e ( t ) E I δ i ( t ) R Queue 2
Queue 1 latency ( σ ) recovery ( γ ) Figure 2:
Queue 1: M t /G/ ∞ Queue from Exposed to Removed.
Queue 2: M t /G/ ∞ Queuefrom Exposed to Infected
Following a previous study of the outset of epidemic (Trapman and Bootsma,2009), we model the arrival rate as a Poisson random variable with a rate pro-portional to the maximum number of encounters between healthy and infectedindividuals. However, the proposed model innovates by considering a time-varying input rate λ ( t ) = βS ( t ) I ( t ) to cover the epidemic’s evolution over time.Observe in Figure 2 that each departure from Queue 1 represents a decrease ofone infected individual. In turn, each arrival corresponds to a new exposition,i.e. a susceptible individual contracting the disease in its latent stage. Thecompartments E , I , and R denote the corresponding SEIR populations.7lso in Figure 2, Queue 2 is a second M t /G/ ∞ queue that covers the latency period σ only. Each departure of this queue increases the infected populationby one.The queues in Figure 2 are drawn from Eick et al. (1993). However, theproposed epidemic model is more general since the input rate λ ( t ) is not adeterministic function of time. Instead, it depends on the number of infectedand susceptible individuals at time t ≥
0. As previously stated, stochastic jumpsin the number of infected individuals depend upon the departure processes ofthe two queues introduced in Figure 2. On the other hand, the susceptiblepopulation is non-increasing over time and decreases with each new arrival.
We start by modelling the arrival process shared by both queues. At anygiven time t ≥
0, the rate of contagion (exposition) is given by λ ( t ) = βS ( t ) I ( t ) , t ≥ , (1)where β is a scalar parameter, S ( t ) and I ( t ) are the sizes of the susceptible andinfected populations at time t , respectively.By design, the output of Queue 2 - in Figure 2 - is the number of newinfections, 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 attime t ≥
0. According to Eick et al. (1993), δ e ( t ) is a Poisson random variablewith mean: δ e ( t ) = E ( λ ( t − σ ) ) . (2)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 thetotal 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: δ i ( t ) = E ( λ ( t − ( σ + γ ) )) . (3)8 .1.2. 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 ≥ t ≥ X ( t ) = ( S ( t ) , E ( t ) , I ( t ) , R ( t ) ) ∈ Ω is the state of process X t , t ≥
0. Bydefinition, X t , t ≥ i ) a new exposition, ii ) a new departure from Queue2 or iii ) a new departure from Queue 1.To describe the dynamics of the process X t , t ≥
0, we make use of the factthat 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 thenext event will be exponentially distributed and the total jump rate at time t ≥ t ) = λ ( t ) + δ e ( t ) + δ i ( t ) . (4)Let { τ , τ , . . . } be the sequence of jumps in the system, with τ ≡ τ k +1 >τ k , ∀ k ≥
0. Now, assume a jump occurs at time t = τ k . Then, the followingholds: P ( X ( t + ) = Y | X ( t ) = ( S ( t ) , E ( t ) , I ( t ) , R ( t ))) , t = τ k )= λ ( t )Λ( t ) if Y = ( S ( t ) − , E ( t ) + 1 , I ( t ) , R ( t )) ) δ e ( t )Λ( t ) if Y = ( S ( t ) , E ( t ) − , I ( t ) + 1 , R ( t )) ) δ i ( t )Λ( t ) if Y = ( S ( t ) , E ( t ) , I ( t ) − , R ( t ) + 1) )0 otherwise . (5)The first expression on the right-hand side of Eq. (5) corresponds to a newcontagion/exposition, which happens at rate λ ( t ), see Eq. (1), and implies thetransference of an individual from the susceptible to the exposed population.The second expression corresponds to a new infection that happens upon thedeparture of an individual from Queue 2, which occurs at rate δ e ( t ), see Eq. (2).9n 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 migratesfrom 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: λ ( t + ) = βS ( t + ) I ( t + ) , (6)where the new values on the right-hand side of the expression above vary asa function of the transition probabilities in (5). Clearly, λ ( t ) remains unal-tered between successive jumps. Given an initial distribution and the transitionprobability in (5), the process X t , t ≥
3. Stochastic stability and reproduction number
Whereas deterministic models rely on the evaluation of trivial equilibriumto derive stability conditions and the so-called reproduction number R (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 condi-tions for a receding epidemic as time elapses. Queuing theory connects thisnotion with the so-called occupation rate, a measure of the input to outputratio as time elapses (Shortle et al., 2018).Consider Queue 1 in Figure 2. Assuming a very large population, the sta-bility of such a queue hinges on the occupation rate (e.g., Shortle et al., 2018): ρ ( t ) = λ ( t ) δ i ( t ) , (7)and can be ascertained if a finite time ¯ t ≥ ρ ( t ) < , ∀ t ≥ ¯ t .This guarantees that the system stabilises and the number of customers in thequeue remains finite. With a finite population, however, stability is guaran-teed because the number of infected individuals will remain within finite, albeitpossibly large, bounds. In that case, we are interested in the trend of Queue10, which indicates whether the epidemic is increasing or receding. Theorem 1below establishes the condition for a receding epidemic. Theorem 1.
Consider the Markov process described in Section 2.1.2 and as-sume that ρ ( t ) < for all t > ¯ t ≥ , with ¯ t < ∞ . Then, it follows that: E (cid:0) E ( t + ) + I ( t + ) | X ( t ) ; t = τ k (cid:1) < E ( t ) + I ( t ) , ∀ t = τ k > ¯ t, k ≥ , (8) where { τ , τ , . . . } is the sequence of jumps in the system, as defined in Section2.1.2. Proof
From (5) we have that: E (cid:0) E ( t + ) + I ( t + ) | X ( t ) ; t = τ k (cid:1) = λ ( t )Λ( t ) ( E ( t ) + I ( t ) + 1) + δ e ( t )Λ( t ) ( E ( t ) + I ( t )) + δ i ( t )Λ( t ) ( E ( t ) + I ( t ) − E ( t ) + I ( t ) + λ ( t ) − δ i ( t )Λ( t ) , for all t = τ k . The last equality holds because Λ( t ) = λ ( t ) + δ e ( t ) + δ i ( t ) - Eq. (4).Now, Eq. (8) follows by assuming ρ ( t ) < t > ¯ t ≥ As a consequence of Theorem 1, we can interpret ρ ( t ) as the reproductionnumber of the system at time t and ρ ( t ) < ∀ t > ¯ t ≥ Theorem 1 provides a basis for developing mitigation strategies to sta-bilise the epidemic, whereas the Markov model in Section 2.1.2 enables us toevaluate the long-term effects of such strategies. Following the literature, weintroduce the control (mitigating action) in the form of non-pharmaceuticalinterventions 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 ) ≤ u ( t ) percentage points in the transmission rate attime t ≥
0, thus resulting in a controlled exposition rate: λ ( t, u ( t )) = β (1 − u ( t )) S ( t ) I ( t ) , t ≥ . (9)Let π = { u ( t ) , t ≥ } be a mitigation strategy and let Π denote the set ofall feasible strategies Π = { π : 0 ≤ u ( t ) ≤ , ∀ t ≥ } . The dynamics of thesystem under any strategy π ∈ Π can be evaluated by running the Markov chain X t , t ≥ λ ( t ) = λ ( t, u ( t )) , t ≥
0. We are particularly interested in the class of stabilising policies Π S ∈ Π,such that: u ( t ) = , if t ≤ ¯ t, min (cid:20)(cid:18) − ¯ ρρ ( t ) (cid:19) , (cid:21) , if t > ¯ t, (10)where ρ ( t ) is the uncontrolled occupation rate in Eq. (7). It is easy to seethat 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 nextSection we will evaluate these policies in the light of the COVID-19 epidemicand explore the stabilising effect of parameters ¯ ρ and ¯ t in Eq. (10). Observe12hat, because the population is finite, the effect of the control is highly sensitiveto and can be limited by a delayed start of mitigation.Observe, however, that the Markov model in Section 2.1.2 is not limited tothe proposed class of mitigation strategies. In fact, one can substitute any ar-bitrary mitigation policy π ∈ Π in Eq. (9) and run the controlled Markov chainto evaluate the effect of such policy. To demonstrate that, our experiments willalso 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 scalelock-down when infections surpass a prescribed upper bound; conversely, miti-gating actions are lifted when infections fall bellow a prescribed lower bound.
4. Numerical Experiments
This Section validates the proposed stochastic SEIR dynamic model in Sec-tion 2.1.2 in the light of data from the COVID-19 epidemic. We replicate theparameters of Tarrataca et al. (2021) concerning the spread of the epidemic inBrazil. The parameters are listed in Table 2.
Table 2: Parameters of the Simulation.
Parameter Description Value P Total population 217 · I (0) Initial number of infected people 2 E (0) Initial number of exposed people 252 S (0) Initial number of susceptible people 217 · − β Transmission rate 2 . · − To complete the necessary data to run the model, we also need the distribu-tions of the latency period σ and of the latency plus infection period ( σ + γ ),which represents the full disease cycle. The distributions are compatible withprevious observations (Backer et al., 2020; Verity et al., 2020) and assumed dis-crete to cope with the typical daily collection of data. Tables 3 and 4 unveil thedistributions of σ and ( σ + γ ), respectively.13 able 3: Latency Period Distribution ( σ ). Day
Probability - 0.0009 0.0056 0.0222 0.0611 0.1222 0.1833 0.2095
Day
Probability · − Table 4: Latency and Infection Period Distribution ( σ + γ ). Day
Prob. - - - - - - - -
Day
Prob. - - 3 . · − . · − . · − Day
16 17 18 19 20 21 22 23
Prob.
Day
24 25 26 27 28 29 30 31
Prob.
Day
32 33 34 35
Prob. . · − . · − . · − . · − The first series of experiments examine policies π = (¯ t, ¯ ρ ) ∈ Π S that followEq. (10) with a fixed target occupation level ¯ ρ = 0 .
95. The objective is toevaluate the effect of delaying the start of mitigating actions by varying thefirst parameter. Cases A-G in Table 5 feature different triggering times forthe mitigating actions, starting with the case where no mitigation is enforced(¯ t = ∞ ). For each experiment, we ran a full realisation of the system for twoyears, starting from the outset of the epidemic. Table 5: Simulated cases.
Case
A B C D E F G¯ t ∞
63 91 98 105 112 126¯ ρ Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 3:
Epidemic evolution for Case A: the unmitigated spread
Figure 3 shows the daily evolution of the epidemic for Case A, i.e. in theabsence of mitigating actions. To facilitate the interpretation, the number ofindividuals in each SEIR compartment is depicted as a percentage of the totalpopulation. One can see that the epidemic spreads rapidly, with the number ofinfected individuals peaking around 64% of the population after 128 days. Theepidemic then starts to shrink due to the decrease in the susceptible population,until it finally subsides around the 165 th day.Observe from Figure 4 that Case B starts the mitigation early in the epi-demic, with low levels of infection. The upper plot conveys full results andin the lower plot we zoom in to detail the evolution of the exposed, infectedand removed populations. The results illustrate the effectiveness of the mitiga-tion policy. Observe that the mitigation policy successfully prevents the spreadof the disease, maintaining 99.6% of the population healthy over the wholetwo-year interval, whereas the removed population reaches 0.4% at the end ofthat interval. Controlling the epidemic, however, demands very high levels ofmitigation that oscillate and stabilise around 0 . Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t)0.0100.0200.0300.0500.0750.1000.2000.3000.400 63 91 200 400 600 730 800 Time (days) popu l a t i on ( % ) exposedinfectedremoved Figure 4:
Epidemic evolution for Case B terministic SEIR-based optimal control approaches also required high levels ofcontrol to curb the COVID-19 epidemic, as can be seen in the relevant work ofPerkins and Espa˜na (2020).In Case C, the mitigation strategy is delayed for an extra three weeks with16
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 5:
Epidemic evolution for Case C respect to Case B. One can see that the mitigating actions are similar, tend tostabilise around the same value and are equally sufficient to curb the epidemic.The difference with respect to Case B is a rather steep increase in the finalnumber of individuals that catch the disease at some point, represented by theremoved population. This number increases from 0 .
4% in Case B to around17% in Case C, thus highlighting the sensitivity of the spread with respect tomitigation delay. The results also show that a three-week delay produces alarge increment in the peak of infections, for even though this peak only reachesaround 1% of the population, it is about four times as large as it would be ifthe mitigation strategy started 21 days earlier.Figures 6 features the results for Case D. Despite an extra delay of one weekwith respect to Case C, it is still possible to flatten the curve and prevent anuncontrollable increase of the epidemic. Nonetheless, the extra week producesa peak of infections five times as large as in Case C, with 5% of people beinginfected 107 days after the outset of the epidemic. The total removed populationskyrockets from about 17% in Case C to about 90% in Case D, thus reinforcing17
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 6:
Epidemic evolution for Case D the importance of early mitigation to prevent the spread of the disease andthe exponential effect of extra delays in the system. Observe that, due to thespread of the disease and the corresponding decrease in the number of susceptibleindividuals, we are able to slowly relax the mitigation until the disease hascontaminated about 50% of the population. As the disease spreads further, weare able to rapidly relax the mitigation whilst still maintaining a steady decreasein the number of infected individuals. The mitigation is completely lifted justbefore the second anniversary of the outset of the epidemic.Figure 7 shows the results for Case E that provide further evidence of thedeleterious effect of mitigation delay. A further delay of one week with respectto Case D results in approximately three times as many infected individualsat the peak, that occurs on day 114 and amounts to around 14% of infectedindividuals. Due to the larger peak, we are able to relax the mitigation earlierand more rapidly, reaching a full relaxation just after seven months. However,since there are still susceptible individuals in the population, mitigation has tore-start from day 396 onwards to avoid a recrudescence of the epidemic. The18
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 7:
Epidemic evolution for Case E final number of removed individuals indicates that the epidemic affects about97.5% of the population within two years.
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 8:
Epidemic evolution for Case F
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 9:
Case G
Finally, Figure 9 presents the results for Case G, which demonstrate thatdelaying the mitigation by 126 days effectively renders it meaningless. Indeed,on day 126, the uncontrolled epidemic is so widespread that the infection levelsare already dropping due to a small number of susceptible individuals remainingin the population. The infection reaches its peak on day 129 when the infection20imultaneously afflicts approximately 64% of the population.
This Section investigates the effect of distinct target occupation rates ¯ ρ in thestabilising policies. To provide an interesting baseline, we start our mitigatingactions on ¯ t = 105, when the epidemic reaches a significant proportion of thepopulation but can still be controlled as in Case E of Figure 7. The set ofexperiments in Table 6 assesses the effect of varying the target occupation levelfrom 60% to 90% in regular intervals. To account for the stochastic fluctuationsand prevent biases, Figures 10 to 13 show the median of 100 realisations of thestochastic system within a two-year interval. Table 6: Second set of experiments.
Case
H I J K¯ t
105 105 105 105¯ δ Figure 10 presents the results for Case H. As one can observe, the numberof infected people peaks around 18% on day 114. As expected, the strongermitigation has no effect in the early stages and is not able to prevent a similarpeak as in the baseline (Case E). However, as time elapses we can notice that theincreased 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. Theincreased control is able to keep the removed population at about 52% at theend of two years, whereas that level was 97.5% in the baseline. It is noteworthythat the mitigation is kept at a high level after the system stabilises to avoidrecrudescence.The results for Case I are shown in Figure 11. The increase in the targetoccupation rate ¯ ρ with respect to Case H has several noticeable effects: a slightdecrease in the mitigating actions; a slower rate of decrease of the infection21 Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 10:
Epidemic evolution for Case H levels and an increase of the removed population. The infection levels decreaseto 3% on day 200 and 72 days later reach 1%. The reduced control levels alsolead to an increase in the removed population, that climbs from 52% in Case Hto about 63% in Case I.The trend of increasing infection levels gathers speed in Case J, as the de-crease in the required mitigating actions leads to about 93% of removals withinthe two-year horizon. On day 200, the number of infections reaches 5% of thepopulation, and 117 days later it reaches 1%. Coupled with the significant de-crease in the susceptible population, which reduces the potential spread, thedecrease in ¯ ρ produces a steep decrease in the control levels. These reach a min-imum of 5% as the susceptible population starts to stabilise. However, increasedlevels of mitigation are required later to avoid a second wave.Finally, Figure 13 presents the results for Case K. As the target occupationlevel approaches that of the baseline (Case E), the overall behaviour becomesquite similar. We can notice a small decrease in the overall removals, that reachabout 3.5% as opposed to 2.5% in the baseline. Despite dropping mitigating22 Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 11:
Epidemic evolution for Case I
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 12:
Epidemic evolution for Case J actions around day 210, the reduced size of the susceptible population causesthe level of infection to fall sharply after reaching 5% on day 237. After the23
Time (days) popu l a t i on ( % ) exposedinfectedremovedsusceptibleu(t) Figure 13:
Epidemic evolution for Case K system stabilises, late mitigation efforts are needed to keep the infection at bayjust before the 500 th day. This Section briefly evaluates the effect of on-off lock-down policies. Pro-posed 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 upperlimit; conversely, all measures are lifted ( u ( t ) = 0 , ¯ ρ = ∞ ) as soon as thesenumbers drop below a lower bound. We simulate the two policies whose upperand lower bounds appear in Table 7. Table 7: Third set of experiments.
Case
L M
Start
Stop Time (days) popu l a t i on ( % ) exposed infected removed susceptible u(t) Figure 14:
Epidemic evolution for Case L
Figure 14 shows that Case L features 12 lock-downs in two years, the firststarting on day 99. Infections recurrently exceed the upper bound in Table 7 andpeak around 4% on day 108; the removed population amounts to 52% at the endof the second year. Consistently with the results reported by Tarrataca et al.(2021), the system alternates between rapid increases and decreases in infectionand the lock-downs become shorter and more spaced over time.A similar behaviour is observed in Figure 15 for Case M. But because theupper and lower bounds are larger, the total number of removals increases from52% to 74% in two years. In addition, the peak of infections observed on day110 increases to 6%. Although Case M has 2 more lock-downs in two years,the total time in lock-down is similar, as illustrated in Figure 16, which showsthe cumulative time in lock-down on the left hand side and the duration ofindividual lock-downs on the opposite half. The mean duration of a lock-downis about 23 days in case L and 21 days in case M, whereas the mean interval25
Time (days) popu l a t i on ( % ) exposed infected removed susceptible u(t) Figure 15:
Epidemic evolution for Case M between lock-downs was 32 days for Case L and 27 days for Case M.
Time (Days) T i m e i n Lo ck do w n ( da ys ) Case LCase M 1820222426 1 2 3 4 5 6 7 8 9 10 11 12 13
Lockdown Sequence D u r a t i on ( da ys ) Case LCase M
Figure 16:
Time in Lock-down in Cases L and M
5. Conclusions
Inspired by the unprecedented COVID-19 epidemic, this paper explored in-novative stochastic modelling alternatives to describe the evolution of an epi-demic by means of the classical SEIR framework. The proposed model innovates26y simultaneously addressing general infection and latency periods, whilst theoriginal methodology ensures that the resulting continuous-time Markov modelremains tractable and useful to support decision making.By drawing a parallel between stochastic stability and the traditional repro-duction number, the paper introduced stabilising strategies that can curb anepidemic under very general conditions, provided that the mitigation is initi-ated in a timely manner. Beyond the academic contributions, we argue that thegenerality of the proposed approach renders it invaluable to support real-worlddecision making in the face of future epidemics. The methodology was validatedusing COVID-19 data from the literature. The results provide a panorama ofinsights into the pros and cons of distinct stabilising mitigation strategies overa two-year horizon.As expected, the experiments illustrate that early intervention is vital toprevent the disease from affecting a large portion of the population. However,to effectively prevent the spread, we need very high levels of mitigation overthe whole two-year horizon. The results also suggest that delayed mitigationleads to a dramatic increase in the overall number of infections. In effect, delaystend to produce persistent infection levels over time, even in the presence ofmitigation, thereby increasing the spread even though the infection curves areeffectively flattened. Beyond flattening the curves, this behaviour suggests thatwe also need to ensure, by acting swiftly, that their summit is tolerable fromboth societal and health care perspectives.The proposed approach leaves many research avenues to be explored in futureworks. One obvious route is to pursue stochastic optimal control policies thatsomehow address the compromise between health care aspects, societal issuesand the economic burden of mitigation strategies. The challenges involve findinga meaningful trade-off among the different elements that decision-makers needto consider, as well as proposing effective formulations that avoid the curse of di-mensionality (Powell, 2011, 2019) to ensure that the problem remains tractable.Another branch goes into developing filtering and analytical approaches to esti-mate a system’s parameters considering that the available information is delayed27nd biased, since dependent on local testing and reporting policies.
Acknowledgements
This study was partly supported by the Brazilian Research Council—CNPq,under grants
References
Allen, L.J.S., 2008. An Introduction to Stochastic Epidemic Models, in: Brauer, F., van denDriessche, P., Wu, J. (Eds.), Mathematical Epidemiology. Springer Berlin Heidelberg, Berlin,Heidelberg, pp. 81–130. doi: .Amador, J., Lopez-Herrero, M., 2018. Cumulative and maximum epidemic sizes for a nonlinearSEIR stochastic model with limited resources. Discrete & Continuous Dynamical Systems - B23, 3137. doi: .Artalejo, J.R., Economou, A., Lopez-Herrero, M.J., 2015. The stochastic SEIR model before ex-tinction: Computational approaches. Applied Mathematics and Computation 265, 1026 – 1043.doi: .Backer, J., Klinkenberg, D., Wallinga, J., 2020. Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China, 20–28 january 2020. Eurosurveillance25. doi: .Barraza, N.R., Pena, G., Moreno, V., 2020. A non-homogeneous Markov early epidemic growthdynamics model. application to the SARS-CoV-2 pandemic. Chaos, Solitons & Fractals 139,110297. doi: .Br´emaud, P., 1999. Gibbs fields, monte carlo simulation, and queues. Springer-Verlag, New York.Britton, T., 2010. Stochastic epidemic models: A survey. Mathematical Biosciences 225, 24 – 35.doi: .Clancy, D., 2014. SIR epidemic models with general infectious period distribution. Statistics &Probability Letters 85, 1 – 5. doi: .Davis, M.H.A., 1993. Markov models and optimization. Chapman and Hall, London.Dike, C.O., Zainuddin, Z.M., Dike, I.J., 2016. Queueing Technique for Ebola Virus Dis-ease Transmission and Control Analysis. Indian Journal of Science and Technology 9.doi: . ick, S.G., Massey, W.A., Whitt, W., 1993. The Physics of the M t /G/ ∞ Queue. OperationsResearch 41, 731–742. doi: .Ferguson, N., Laydon, D., Nedjati Gilani, G., Imai, N., Ainslie, K., Baguelin, M., Bhatia, S.,Boonyasiri, A., Cucunuba Perez, Z., Cuomo-Dannenburg, G., Dighe, A., Dorigatti, I., Fu, H.,Gaythorpe, K., Green, W., Hamlet, A., Hinsley, W., Okell, L., Van Elsland, S., Thompson, H.,Verity, R., Volz, E., Wang, H., Wang, Y., Walker, P., Winskill, P., Whittaker, C., Donnelly,C., Riley, S., Ghani, A., 2020. Report 9: Impact of non-pharmaceutical interventions (NPIs) toreduce COVID-19 mortality and healthcare demand. Technical Report. Imperial College London.doi: .G´omez-Corral, A., L´opez-Garc´ıa, M., 2017. On SIR epidemic models with generally distributedinfectious periods: Number of secondary cases and probability of infection. International Journalof Biomathematics 10, 1750024. doi: .Hethcote, H.W., 2000. The Mathematics of Infectious Diseases. SIAM Review 42, 599–653.doi: .Kantner, M., Koprucki, T., 2020. Beyond just “flattening the curve”: Optimal control of epi-demics with purely non-pharmaceutical interventions. Journal of Mathematics in Industry 10,23. doi: .Kermack, W.O., McKendrick, A.G., Walker, G.T., 1927. A contribution to the mathematical theoryof epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of aMathematical and Physical Character 115, 700–721. doi: .Lef`evre, C., Simon, M., 2020. SIR-Type Epidemic Models as Block-StructuredMarkov Processes. Methodology and Computing in Applied Probability 22, 433–453.doi: .Lopez-Herrero, M., 2017. Epidemic Transmission on SEIR Stochastic Models with Nonlinear Inci-dence Rate. Mathematical Methods in the Applied Sciences 40, 2532–2541. doi: .L´opez-Garc´ıa, M., 2016. Stochastic descriptors in an SIR epidemic model for heterogeneous individ-uals in small networks. Mathematical Biosciences 271, 42 – 61. doi: .Meyn, S.P., Tweedie, R.L., 1993. Markov Chains and Stochastic Stability. Springer-Verlag, NewYork.Perkins, T.A., Espa˜na, G., 2020. Optimal Control of the COVID-19 Pandemicwith Non-pharmaceutical Interventions. Bulletin of Mathematical Biology 82, 118.doi: .Powell, W., 2011. Approximate Dynamic Programming Solving the Curses of Dimensionality. JohnWiley & Sons, Inc., New Jersey, USA.Powell, W., 2019. A Unified Framework for Stochastic Optimization. European Journal of Opera-tional Research 275, 795–821. doi: . oss, R., 1916. An application of the theory of probabilities to the study of a priori pathometry-partI. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematicaland Physical Character 92, 204–230. doi: .Shortle, J., Thompson, J., Gross, D., Harris, C., 2018. Fundamentals of Queueing Theory. WileySeries in Probability and Statistics. 5 ed., Wiley, New York. doi: .Tarrataca, L., Dias, C.M., Haddad, D., Arruda, E.F., 2021. Flattening the curves: on-off lock-downstrategies for COVID-19 with an application to Brazil. Journal of Mathematics in Industry 11,2. doi: .Trapman, P., Bootsma, M.C.J., 2009. A useful relationship between epidemiology and queueingtheory: The distribution of the number of infectives at the moment of the first detection. Math-ematical Biosciences 219, 15 – 22. doi: .van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemicequilibria for compartmental models of disease transmission. Mathematical Biosciences 180, 29– 48. doi: https://doi.org/10.1016/S0025-5564(02)00108-6 .Verity, R., Okell, L.C., Dorigatti, I., Winskill, P., Whittaker, C., Imai, N., Cuomo-Dannenburg,G., Thompson, H., Walker, P.G.T., Fu, H., Dighe, A., Griffin, J.T., Baguelin, M., Bhatia,S., Boonyasiri, A., Cori, A., Cucunub´a, Z., FitzJohn, R., Gaythorpe, K., Green, W., Ham-let, A., Hinsley, W., Laydon, D., Nedjati-Gilani, G., Riley, S., van Elsland, S., Volz, E., Wang,H., Wang, Y., Xi, X., Donnelly, C.A., Ghani, A.C., Ferguson, N.M., 2020. Estimates of theseverity of coronavirus disease 2019: a model-based analysis. The Lancet Infectious Diseasesdoi: .˙I¸slier, Z.G., G¨ull¨u, R., H¨ormann, W., 2020. An exact and implementable computation of the finaloutbreak size distribution under Erlang distributed infectious period. Mathematical Biosciences325, 108363. doi: ..