key: cord-1034109-ef1s7zzi
authors: Kalbaugh, David V.
title: Probabilistic Model for Control of an Epidemic by Isolation and Quarantine
date: 2021-04-23
journal: Bull Math Biol
DOI: 10.1007/s11538-021-00897-1
sha: bb4b6ca9b2c1a38fa2059f21b2e13e197ecb2d11
doc_id: 1034109
cord_uid: ef1s7zzi

Assuming a homogeneous population, we apply the mass action law for rate of new infections and a second-order gamma distribution for removal probability to model spread of an epidemic. In numerical examinations of higher-order gamma distributions for removal probability, we discover an unexpected pattern in maximum fraction of population infected. We develop from first principles of probability an eighth-order system of ordinary differential equations to model effects of isolation and quarantine. We derive analytical expressions for reproduction numbers modeling isolation and quarantine when applied separately and together and verify them numerically. We quantify strength and speed required of these interventions to contain epidemics of varying severity and examine how their effectiveness depends on when they begin. We find that effectiveness is highly sensitive to small changes of intervention strength in a critical region. Finally, adding two more differential equations to capture natural population dynamics, we calculate endemic disease equilibria when affected by isolation and examine dynamics of coming to an equilibrium state.

Mathematical epidemiology has a long history. In 1760, Daniel Bernoulli analyzed inoculation against smallpox, where he developed what is usually described as the first model in the subject (Bernoulli 1766; Blower 2004; Heesterbeck and Roberts 2015) .

On developments directly relevant to the matter at hand, in 1906 Hamer proposed a mass action law, in which rate of new infections in a disease outbreak is proportional to the product of the number of susceptible individuals and the number of infected individuals. This model has been widely used since that time (Hamer 1906; Brauer 2017) . Kermack and McKendrick (1927) introduced compartments to model spread of infectious diseases. A special case of their general model is referred to as the SIR model, with compartments S (susceptible but well), I (infected and able to transmit the disease) and R (removed by immunity, recovery, death or quarantine). In this model, people who recover are permanently immune. The SIR model incorporates the mass action law and has become the starting point for a wide variety of compartmental models to the present day. The SIR model is deterministic, meaning it incorporates no element of chance. We will discuss the SIR model more, later in this article.

Beginning in 1928 and continuing through the 1930s and 1940s, L. J. Reed and W. H. Frost lectured at what is now the Johns Hopkins Bloomberg School of Public Health and developed the stochastic Reed-Frost model, which has been widely used and from which many extensions have been formulated (Lessler and Cummings 2016; Brauer 2017) . Their work was published by their students (Abbey 1952; De Maia 1952) . Stochastic models include the element of chance. We will discuss the Reed-Frost model more, later in this article.

Isolation is removal of an infected person from public interaction. Quarantine is removal of a segment of the population, regardless of health condition, usually by government order. A deterministic study by Hethcote et al. (2002) modeled compartments SIQR and SIQS. In the latter, people who recover are not immune but rather return to susceptibility. The compartment labeled Q was for isolation, in that transitions to Q were from I only. The analysis found conditions under which disease can be endemic, i.e., have a permanent presence. In some cases, the presence was a stable, constant equilibrium. In others, it was eventually periodic. Gumel et al. (2004) authored a study on strategies for controlling SARS (severe acute respiratory syndrome) outbreaks. They used a deterministic SEQIJR model, where Q and J were for quarantine and isolation and E labeled a compartment for the exposed but asymptomatic. In the Hethcote et al. study above, the Gumel et al. study and most compartment-based studies, transitions between compartments are modeled as first-order, linear, single-parameter rates except for rate of new infections, which is modeled by various versions of the mass action law. In the Gumel et al. study, transitions to the quarantine compartment were from the exposed compartment only, which implies they were based on contact tracing. The study compared model predictions to actual outcomes in a number of geographic areas and found good agreement. It also showed the benefit of fast response in both isolation and quarantine together.

The literature has many articles analyzing isolation and quarantine that are based on compartmental models. Many are deterministic (e.g., Castillo-Chavez et al. 2003; Sahu and Dhar 2015; Denes and Gumel 2019) and some are stochastic by virtue of adding Wiener processes (e.g., Li et al. 2018) . The latter require multiple solutions of the equations to produce multiple realizations of an epidemic in order to yield reliable "expected value" results (Monte Carlo simulation).

All the isolation and quarantine models discussed above assume a homogeneous population. Other studies using the same assumption but not employing compartmental models include Day et al. (2006) , which used a discrete time branching process, and Fraser et al. (2004) , which utilized partial differential equations of the transport type.

Even if populations were homogeneous, compartmental models are inadequate for the very early dynamics of epidemics (Diekmann et al. 2013; Brauer 2017) . Important fluctuations in initial spread are represented much better by network branching processes.

Spatial, age and demographic heterogeneity are real-life conditions that require additional model sophistication. The professional standard is set by immense computer programs that account for actions of millions of people individually: their age, demographic status and location, with real data on such matters as airline schedules from city to city. The programs are stochastic but so large that only a few realizations can feasibly be produced, even on supercomputers. Halloran et al. (2008 ), Ferguson et al. (2006 ), German et al. (2006 and Lewis et al. (2007) are examples.

In today's world, the sophisticated technologies described above are essential. But simple models have a place. With minimal time and effort, simple models can enable preliminary assessment of a wide range of strategic policy options and point the way to more in-depth studies. They can be useful introductory educational tools for those entering the epidemiology profession. Furthermore, even for basic subjects such as isolation and quarantine, a simple model can quantify, with a broad brush, effort required to contain an epidemic and clearly demonstrate how deleterious delay and noncompliance can be. Finally, it can identify endemic equilibria, conditions portending a long-term struggle against the disease.

This article develops and employs such a model. Section 2 develops the model from first principles of probability as it applies to unmitigated spread of an epidemic. Section 3 compares that model to the SIR equations and finds an interesting ratio in maximum fraction of population infected, which Sect. 4 explores further. Applying first principles of probability, Sect. 5 extends the model to include effects of isolation and quarantine. Section 6 employs the model to (a) determine requirements on isolation and quarantine to contain an epidemic, (b) quantify mitigation effectiveness as it depends on when intervention begins and (c) measure consequences of efforts falling short. It finds that intervention effectiveness is highly sensitive to small changes of effort in a critical region. Section 7 extends the model once more to examine endemic disease equilibria when affected by isolation. Finally Sect. 8 reviews main points of the article and summarizes advantages and disadvantages we perceive in this model with respect to others in the literature.

Letŝ j (t) = conditional probability that individual j in a population of N people is susceptible to the disease but has not contracted it after time t since its outbreak, given the sequence of contacts and times t k , k = 1, … ,n j (t). i j (t k ) = probability that the person with whom individual j interacted in his or her k th interaction since disease outbreak was infected at time of interaction. c = probability that the disease is transmitted in a single interaction between a susceptible-but-well person and an infected one. n j (t) = number of interactions individual j has had with other people by time t.

Then, ŝ j (t) is given by the product of the probabilities of the individual not being infected on each of his or her n j (t) interactions:

We are assuming for simplicity that no one is initially immune. Equation (1) is a form of the Reed-Frost model. The variables n j (t) , ŝ j (t) and ̂i j (t) are stochastic processes. The number n j (t) takes on only integer values, increasing by one at randomly occurring interactions. Any realization of ŝ j (t) is constant between interactions, where it takes discontinuous steps down. From this, we form a deterministic, differentiable function by taking expected values.

Averaging over all possible interactions, the individuals involved and when,

We assume that interaction times follow a homogeneous Poisson process with growth rate a , assumed constant. Let X =n j (t + Δt) −n j (t) . Then X is Poisson distributed with mean aΔt . We also assume that i j (t) changes little in the small Δt interval. Then we can approximate

(1 − ci j (t)) k (aΔt) k k! = exp(−ci j (t)aΔt).

Assuming homogeneity of the population, this becomes Then Taking the limit as Δt → 0, This is the aforementioned mass action law. We have developed a well-known deterministic model from a well-known stochastic one. The fundamentals of this relationship are discussed in basic texts (e.g., Diekmann et al. 2013) . We have included a bare sketch of the derivation to show the probabilistic underpinnings of our model. For a detailed rigorous proof, see Kurtz (1981) . The product ac is often denoted , and we will follow that custom henceforth.

Next, we must derive an equation for i(t) . Let r(t) be the probability a given person has been infected at some time in the epidemic but has been removed via recovery or death by time t . We assume that a person who recovers is permanently immune. To compute r(t) , we begin with the probability distribution which has probability density 2 R (t − ) exp(− R (t − )) . This says an infected person is more likely to be removed at a time 1∕ R after onset of infection than at any other time. Equation (5) is the second member of the Erlang family of distributions, which is the subset of the gamma family having integer order.

Let p(t) be the probability a given person was infected at some time prior to t in the epidemic. Then

The state equation for p is Also, s j (t + Δt) = s j (t) exp(−ci j (t)aΔt).

(3) s(t + Δt) = s(t) exp(−ci(t)aΔt).

(4) ds dt = −acsi. (6) s = 1 − p.

(7) dp dt = si.

(8) i = p − r if we assume that deaths due to the epidemic are not an appreciable fraction of the population. It is convenient to say here only that a nonconstant circulating population is a complication and to defer further explanation of the need for this assumption to Sect. 5. From Eq. (5) and definition of p:

Integrating by parts:

Numerical values for r(t) as given by Eq. (10) can be computed most readily by the following general process. Anticipating needs later in this paper, we consider an equation of the form Differentiating:

Denote the second term on the right-hand side of Eq. (12) as 1 . Then, Eq. (12) becomes Differentiating 1 :

Denote the second term on the right-hand side of Eq. (14) as 2 . Continuing in this way

Numerically integrating the simultaneous ordinary differential Eqs. (13), (15), (16) and (17) gives numerical values to 0 (t) as given by Eq. (11).

Translating Eq. (10) into the general form in Eq. (11) results in the following equations for r(t):

Equations (6), (7), (8), (18) and (19) constitute our basic model, applying to spread of an unmitigated epidemic. For simplicity and consistency, we will refer to it as the basic probabilistic model, although at this point, as we will see in the next section, it differs from the SIR model only in the equation for removal probability. The name probabilistic model will be more fitting when, using first principles of probability, we extend it to include isolation and quarantine in Sect. 5 and natural population dynamics in Sect. 7.

The SIR equations can be written (See, e.g., Hethcote 2000.) It is not difficult to show that the SIR model assumes and for brevity we omit the proof. Equation (23) is the first member in the Erlang family of probability distributions. The probability density for the distribution in Eq. (23) is exp(− (t − )) , which says that an infected person is more likely to recover immediately after infection than at any later time. This is unrealistic, as

discussed, e.g., in Diekmann et al. (2013) , so we prefer to use the probability distribution in Eq. (5).

The SIR equations are attractive for their simplicity and utility. Two important measures of epidemic severity are (a) s ∞ , probability that a person chosen at random from the population is susceptible but well at the end of the epidemic and (b) i MAX , maximum probability of that person being infected (largest fraction of the population infected at any one time) during the epidemic. Derivations of the SIR model's well-known equations for these measures are included here for convenience.

Setting di∕dt = 0 in Eq. (21), we find that s = ∕ is the fraction of people susceptible but well when the number of people infected in the epidemic reaches its peak. Using Eqs. (20) and (21), we find that Integrating:

Equation (26) applies beyond the SIR model, as we will discuss. To compare the SIR and the basic probabilistic model quantitatively, we must find common ground, setting each model's parameters so they represent the same severity of epidemic. Common ground is found in the basic reproduction number R 0 . R 0 can be defined as the expected number of cases directly generated by one infected person in a population where all other individuals are susceptible to infection. Mathematically, For the SIR model, as is well-known and easy to prove:

For the basic probabilistic model:

Prob{infected person is in circulation at time t} ⋅ average contact rate at time t ⋅ Prob{transmitting infection on one interaction at time t}dt

So for the basic probabilistic model, Therefore, to compare the SIR and basic probabilistic models we set R = 2 . Figure 1 graphs s ∞ and Fig. 2 graphs i MAX for the two models as a function of basic reproduction number. Numerical integration was accomplished with the Adams-Moulton algorithm.

There has been interest in the generality of the final size formula, which is Eq. (26) with ∕ replaced by R 0 . (Final size is defined as 1 − s ∞ .) For example, Ma and Earn (2006) demonstrate that the formula applies to three SIR variations. Figure 1 shows that the final size formula applies as well to the basic probabilistic model. In fact, Anderson and Watson (1980) proved that, for SIR-like models, removal probabilities given by gamma distributions of any order yield the same final size. Figure 2 shows that i MAX for the probabilistic model can be considerably larger than that for the SIR model. Interestingly, to the accuracy of our numerical integration, the ratio of i MAX for the second-order removal distribution is 4/3 that for the SIR model as R 0 approaches one. This curious result will divert us temporarily to a mathematical sidebar in the next section. The ratio is not greatly different from 4/3 when R 0 is much larger; e.g., the ratio is 1.300 when R 0 is 2.75.

We consider the SIR model with higher-order members of the Erlang family of probability distributions for the removal process. That is, we consider equation sets of the form which produce a SIR-like model with removal probability density where t − is time since onset of infection and 1∕ the average time from infection to removal. We have computed results for n up to and including 5 and found that within the accuracy of our numerical integration:

The limit of the ratio of i MAX produced by the gamma distribution of order n to that given by Eq. (25) is 2n∕(n + 1) as ∕ goes to one.

We conjecture that this is mathematically true for all integer orders.

Our objective here is to develop a probabilistic model of quarantine and isolation in an epidemic. Section 6 will employ the model to quantify important aspects of their effectiveness.

We are modeling a homogeneous population but ought to take into account that governments differ and orders to quarantine can come from anywhere within the hierarchy of nation, state, county and city. We can do this by assuming that the order comes with a probability distribution. Compliance to the order will also be distributed and we combine the two events into one process modeled by.

where q A0 is the probability the person will eventually be quarantined and 1∕ A the most probable time after epidemic outbreak that he or she enters quarantine. We believe that modeling both strength and speed of the process is framework for more realistic representation.

Let q AS (t) denote the probability a randomly drawn person is susceptible but well and quarantined. Using Eq. (6), we can compute q AS (t) by

We implement Eq. (32) in the ODEs, using Eq. (31), by

Prob{person drawn at random from within the population is in quarantine at time t after start of the epidemic} =

Equation (33) assumes a susceptible-but-well person once quarantined remains so for the duration of the epidemic.

Let s(t) denote the probability a person drawn at random from the population at time t from start of the epidemic is susceptible but well and not quarantined. From Eqs. (6) and (32) Integrating Eq. (32) by parts, Eq. (34) can be rewritten as which will be useful later on.

In addition to responding to an order, a person may also be removed by the process of isolation after becoming infected at time and recognizing that fact through testing or onset of symptoms. Denote by q B (t| ) the probability of removal by time t via this second process, given infection at time and that the person has not been removed by recovery, death or quarantine. We model q B (t| ) as:

where q B0 is the (conditional) probability of a person eventually going into isolation on learning of his or her infection and 1∕ B is the (conditionally) most probable time the person enters isolation after the infection has occurred. Introducing this process will later cause us to amend the definition of q A (t).

Let i(t) denote the probability that a person drawn at random from the population has been infected but at time t has not been removed by recovery, death, quarantine or isolation. Putting quarantine aside for a moment and recalling Eq. (5), Integrating by parts Substituting Eqs. (5) and (36) into the above and doing a little work lead to

Translating Eq. (39) into the general form in Eq. (11) results in the following equations for i B (t):

We will need to determine the probability q B (t) a person drawn at random from the population is in isolation at time t . We assume that only a person who is currently infectious is subject to isolation. A person goes into isolation based on testing and/or obvious symptoms. Testing and/or disappearance of symptoms should also prove recovery. It is undesirable from both public health and individual liberty points of view for a recovered (and by assumption therefore immune) person to be isolated. Hence, we can describe q B (t) as where r(t| ) is given by Eq. (5), q B (t| ) by Eq. (36), and we will reconsider ṗ( ) shortly. Equation (41) can be rewritten as Integrating by parts, we find that the probability of being put into isolation can be written where i B is given by Eqs. (40) and r by Eqs. (18) and (19). Note that substituting Eq. (42) into Eq. (38) leads to the simpler and more intuitive Equation (43) applies when isolation is implemented alone. When quarantine and isolation are applied together, denoting the probability of a randomly drawn person being infected and in quarantine as q AI : (46) informs us that we must amend the definition of q A (t) alongside Eq. (31). It is not the probability of a given person being in quarantine at time t ; it is the conditional probability given that the person has not previously been isolated.

What effect do quarantine and isolation have on the equation for rate of new infections, i.e., on the equation for p(t) given by Eq. (7) when quarantine and isolation are absent? First, it is clear we must replace i with i∕(1 − q) where q is the probability a person drawn at random from the population is in either isolation or quarantine. (We have not written an equation for q here and will not do so. It would require three more ODEs and in the end we will not need to calculate it.) To motivate the replacement, consider again the person-to-person interactions Sect. 2 analyzed. There, i denoted the probability a person one interacts with is infected, but it also represented the probability a person drawn from the population as a whole is infected. In Sect. 2, s + i + r = 1 . Here, in Sect. 5, i and s include the condition that the person is not in isolation or quarantine. What is important is the probability of interacting with an infected person within the circulating population. Our analysis has not required us to determine the probability a person was removed and not in isolation or quarantine, but if it had, and calling that probability r NQ , we would find s + i + r NQ = 1 − q . To make the probabilities of compartments of people one interacts with sum to unity, we must divide them by 1 − q . Hence, the probability that a person one interacts with is infected is i∕ (1 − q) . Hethcote et al. (2002) and Li et al. (2018) applied this also in SIR-like analyses, with the premise of a constant contact rate, and called si∕(1 − q) quarantine-adjusted incidence. Similar reasoning explains the assumption in Sect. 2 that deaths due to the epidemic are not an appreciable fraction of the population. Without the assumption, Eq. (8) is inconsistent with Eq. (1).

It is also clear that no adjustment need be made to s . The equation for rate of new infections needs the probability a person is susceptible but well and not in isolation or quarantine.

Does the average rate of contacts a , and therefore , require an adjustment? Here the answer is not so clear. We believe there is basis for decision in neither probability theory nor biology of disease and we must make a judgment based on common sense, which is not easy. What, e.g., is the effect of social distancing orders? Hethcote et al. (2002) considered two alternatives. The first, as mentioned above, was to leave contact rate unchanged from its value in an unmitigated epidemic. The second assumed contact rate is reduced in proportion to density of circulating population. We believe that to analyze both here would be unwieldy. Either choice could be defended. Leaving contact rate unchanged leads to simpler equations and provides bounds on performance. However, it leads to equations that we believe undervalue quarantine. Assuming contact rate is reduced in proportion to density of circulating population ascribes, we believe, more appropriate value to quarantine. We reduce contact rate a to a(1 − q) , which modifies to

(1 − q).

Equations (18) To complete the probabilistic model for quarantine and isolation, we derive formulas for reproduction numbers applicable to isolation alone, quarantine alone, and isolation and quarantine together. These expressions are an efficient means of determining minimum effort needed to contain an epidemic, although we will find that containment is achieved only if the processes are begun early enough.

For isolation alone, we return to the definition: expected number of cases directly generated by one infected person in a population where all others are susceptible to infection. The infected individual's interaction with others is unaffected but his or her time in circulation is reduced. Given that the individual was

infected at time = 0 , the probability he or she is still infectious and circulating at time t is (1 − q B (t|0))(1− r(t|0 )) . From Eqs. (5), (27) and (36), We arrive at where Recalling we find Now R 0 = 2 ∕ R is the basic reproduction number for the unmitigated epidemic. Hence, Eq. (52) can be written In the next section, ODE results show that setting R B = 1 in Eq. (54) is remarkably consistent in predicting the strength ( q B0 ) and speed ( B ) required to contain an epidemic, provided isolation is begun early enough. Section 7 further confirms its validity.

Quarantine is much less likely to be applied alone than is isolation, because of its impact on the public, but we consider the possibility here nevertheless.

It is worthwhile to let quarantine reproduction number depend on when quarantine is initiated. In the following discussion, we will be anticipating some of the numerical results of the next section. If quarantine starts early enough, the process can essentially be completed by the time the fraction of infected people becomes appreciable and, using Eq. (35), we can write

But p << 1 , so Similarly, using Eq. (47), with q B = 0, Then Eq. (48) becomes which is just the constant (1 − q A0 ) times the ṗ equation for initial moments of the unmitigated epidemic. Hence, the effect of quarantine with an early start is to reduce to (1 − q A0 ) . There being no dynamic remaining other than infection and recovery, we know that an epidemic will occur in this case if (1 − q A0 )R 0 > 1 . Hence, denoting the reproduction number of quarantine started early as R AE :

The next section will show that this simple equation is useful when quarantine starts early enough.

If quarantine starts late, and we will quantify "early" and "late" in the next section, dynamics of quarantine are simultaneous with those of infection and removal. In this case, we turn once more to the definition of reproduction number based on one infected individual in a population of susceptible but well people. Using Eq. (27), we find A straightforward calculation leads to

When isolation and quarantine are combined, we again have different formulas for early and late start of quarantine. In an early start, A and B processes are uncoupled in time and it is easy to see that the reproduction number for isolation and quarantine combined is For the late start, given the forms of R AL and R B above we might expect the combined reproduction number to be in the form where t B is given by Eq. (53), u AL and v AL by Eqs. (61), and w ABL and x ABL are to be determined. A tedious calculation starting with confirms the form in Eq. (63) and determines that For convenience of the reader, we collect and display below the probabilistic model for control of an epidemic by isolation and quarantine, with exception of the long equations for R ABL above.

(61) (Eqs. 63 and 65) 

We gauge the effectiveness of control efforts by measuring the improvement in s T∞ , the fraction of population not infected during the epidemic, and i TMAX , the largest fraction of population infected at any one time. We conduct analysis at four levels of epidemic severity, for which the probabilistic model predicts s ∞ and i MAX as shown in Table 1 when the epidemic is unmitigated. We base the four epidemic cases on specific values of s ∞ . The unit of time throughout is 2∕ R , the average time of removal in an unmitigated epidemic. For numerical integration, we employ the Adams-Moulton algorithm. Figure 3 describes minimum requirements on q B0 and B to contain an epidemic with isolation alone, on condition that the process is initiated sufficiently early. Figure 3 is based on setting the reproduction number R B in Eq. (56) equal to unity. When values for strength ( q B0 ) and speed ( B ) from Fig. 3 were input to the ODEs and initial condition set at p(0) = 10 −7 the resulting s T∞ values were always between 0.99953 and 0.99959. When p(0) = 10 −3 , those same q B0 and B values produced s T∞ 's varying between 0.95535 and 0.95576, again remarkably consistent but a low bar for containment.

All the intervention reproduction numbers we consider here have this property: If we (a) find strength and speed combinations that yield the reproduction number equal to one and (b) input those strength and speed values into the ODEs, the resulting s T∞ values cluster closely about one point. For some p(0) that cluster point is close enough to unity that the epidemic can be considered contained. As p(0) increases, the cluster point decreases.

We can measure the dependence of the cluster point, call it degree of mitigation ( d M ), on p(0) . Again, in short, d M is the average s T∞ achieved by intervention effort that yields its reproduction number equal to one. Figure 4 graphs d M as a function of p(0) for each of the intervention reproduction numbers considered here. The s T∞ ceiling, 1 − p(0) , is also plotted.

The particular numbers one calculates to create Fig. 4 depend to a degree on parameter ranges one considers: on basic reproduction number R 0 , speeds A and B , and initial condition p(0) . In this study, we arbitrarily selected 1.1886 ≤ R 0 ≤ 2.5582 Table 1 ), 1 ≤ A ≤ 4 , 1 ≤ B ≤ 4 and 10 −7 ≤ p(0) ≤ 10 −3 . Within this set, for each intervention reproduction number R I in Fig. 4 except R B , we can state the following: For each R 0 and each specific effort (q A0 , q B0 , A , B ) satisfying R I = 1 , there exists a z > 1 such that, with these inputs, the ODEs produce and z is nearly independent of p(0). Hence, for each R I but R B , Table 2 lists z for each R I other than R B . Furthermore, our numerical analysis shows that, for R B , where, if expressed as (p(0)) k , k = 1∕2 to better than one part in a thousand.

Returning to discussion of requirements, it is simpler to state than to graph requirements to contain an epidemic with quarantine alone when the process starts early. Setting R AE = 1 in Eq. (58), the minimum q A0 to contain an epidemic is then R AE says quarantine performance is invariant to A , while R AL compares A to removal rate R . The growing dependence on A is gradual. What is an "early" quarantine start? Equation (68) means that both R AE (Eq. 58) and R AL (Eqs. 60 and 61) are consistent throughout the stated p(0) range. The key factor in deciding which to use is how much their s T∞ scatter about their d M . For each intervention reproduction number R I , call

where extrema are again taken over the above-defined set at a particular p(0) . We find that R AL is much more accurate everywhere, but we would prefer, if possible, to use the much simpler R AE . Scatter in s T∞ on the order of 5 ⋅ 10 −6 ( p(0) = 10 −7 ) is, when viewed in the absolute, not a problem for practical purposes. Scatter on the order of 0.05 ( p(0) = 10 −3 ) certainly is. Acknowledging its subjective elements, we offer our opinion that "early" quarantine start ends and "late" start begins at about p(0) = 10 −5 . Figure 5 graphs minimum requirements on q A0 and A to contain an epidemic with quarantine alone, based on setting the "late start" reproduction number R AL in Eq. (60) equal to unity. Comparing Fig. 5 to Fig. 3 , we see that, to contain an epidemic of a given severity, quarantine requires less strength and less speed than does isolation. (Of course, quarantine has a much greater impact on the public.)

When isolation and quarantine are combined, a pertinent issue is allocation of effort between the two. How is burden of performance to be divided? We will lay out examples of choices.

Simplicity of the reproduction number formula for isolation combined with quarantine started early makes allocation for this case an easy task. Allocation is more complicated for late quarantine start, and it is perhaps a more realistic scenario, so we devote our attention to this case.

Figures 6 and 7 quantify strength and speed of the two processes needed to contain an epidemic, on condition that the processes start early enough. For brevity, we show only results for cases II and III. Figures 6 and 7 are based on setting R ABL = 1 , using Eqs. (63) and (65).

It goes without saying that on-the-ground realities drive what can be done; practicalities of implementation are primary. For instance, are means in place for testing and contact tracing necessary for isolation? What is the public's commitment to quarantine? Figures 6 and 7 provide broad-brush quantitative insight into potential of the combined measures. Figure 8 plots s T∞ as a function of q 0 with 0 as parameter for each of the epidemic cases considered, at p(0) = 10 −4 . Figure 8 shows that mitigation achieved by strength q 0 just above a critical value q 0C is much greater, though short of containment, than that achieved by strength just below it. The critical strength is roughly two-thirds of the containment requirement derived from R ABL = 1 , which we denote as q 0R . Table 3 compares q 0C to q 0R for the four epidemic cases, for 0 = 4 and p(0) = 10 −4 .

In Fig. 8 are bold circles marking points where a curve's abscissa meets a containment requirement q 0R . Notice that the circles cluster at about s T∞ = 0.9995. These are but a few of the points that went into calculation of d M for R ABL .

For brevity, we omit graphs of quarantine alone and isolation alone in the format of Fig. 8 . We can summarize important aspects of them: The graph for quarantine alone shows discontinuities similar to those in Fig. 8 . The graph for isolation alone, on the other hand, has smooth, gently sloping curves. The graph of isolation alone also shows that, while d M is only 0.98590 at p(0) = 10 −4 , plausibly achievable increases in q B0 from q B0R values can raise s T∞ to much higher levels for epidemic cases I and II. In case I, e.g., with B = 4, increases in q B0 from the q B0R value of 0.2678 to 0.320 and 0.420 raise s T∞ to 0.997 and 0.999, respectively. Creating Fig. 8 is probably the best way to find levels of intervention effort that achieve desired s T∞ values when d M is too low, i.e., when intervention has begun too late.

In Fig. 8 , s T∞ is highly sensitive to small changes in q 0 . If s T∞ is a continuous function of q 0 , then for every > 0 there exists a > 0 such that if | | Δq 0 | | < then | | Δs T∞ | | < . We conducted numerical experiments at the critical q 0 points to quantify this. An example result: for Case III and 0 = 4 , as | | Δq 0 | | decreased from 10 −3 to 10 −6 Δs T∞ decreased only from 0.24620 to 0.23521. But despite the extremely slow convergence, we believe the function is continuous and, in fact, differentiable. As mentioned earlier, we used the Adams-Moulton (A-M) numerical integration algorithm for this study. To verify that the high sensitivity of s T∞ to changes in q 0 we see in Fig. 8 is not an artifact of the algorithm, we also used the Runge-Kutta (R-K) and backward differentiation formula (BDF) algorithms to integrate the ODEs and recreate Fig. 8 . We found the same high sensitivity with all three algorithms. Furthermore, the difference between q 0C values computed by the A-M algorithm and R-K algorithm was less than 6 ⋅ 10 −6 , and difference between q 0C values computed by A-M algorithm and BDF algorithm was less than 1.3 ⋅ 10 −5 . Differences in computed Δs T∞ were also minimal. Figure 9 plots i TMAX as a function of q 0 with 0 as parameter for epidemic cases III and IV. Discontinuities are again evident, at the same critical values of q 0 . The figure shows that half measures leave the public facing a still-serious epidemic.

Our final objective is to examine endemic equilibrium states. We will calculate equilibrium values for the important state variables and examine dynamics of a disease coming to equilibrium. We assume that isolation, but not quarantine, is in effect. Maintaining quarantine over a significant portion of the population for a very extended period is unlikely.

For a model to exhibit an endemic state, it must include the population's natural dynamics: births, deaths and net immigration. It then becomes necessary for us to use, instead of probabilities, state variables S(t) , the expected number of people in the population susceptible but well at time t ; I(t) , the expected number of people infected and not in isolation at time t ; and R(t) , the expected number of people in the population who have been infected but have recovered at time t . We continue to assume that the number of people who die from the disease is not an appreciable fraction of the population. Denote the number of births and net immigrants per unit time, assumed constant, as B . Let the probability a person drawn at random from the population is dead of natural causes by time t given he or she is alive at time be Let the rate of new infections be given by 

(79) where R B is the isolation reproduction number given by Eq. (54). Note that the threshold for an endemic state in an epidemic with isolation is R B > 1 . Although the probabilistic model has the added complexities of a second-order removal

.

probability and an isolation process, it has led to the well-known SIR equations for the equilibrium state, with R B replacing R 0 . (See, e.g., Diekmann et al. 2013.) How does a disease become endemic? Only under special conditions, it appears. The disease must first have been epidemic, to drive the number of susceptible people below S * . The equilibrium position is at the threshold for herd immunity, when the portion of the population immune to the disease suffices to negate an epidemic. (See, e.g., Hethcote 2000.) If an incursion of infected people occurs when S > S * and mitigation is inadequate, an epidemic inevitably follows, driving S past S * to the S ∞ concomitant with the starting point and I to zero. Then S begins a disease-free course to the natural capacity B∕ D .

If the disease incursion occurs when S is well below S * , the outbreak is relatively quickly extinguished and S continues its path to B∕ D . Only if the incursion happens as S nears S * from below, as in Fig. 10 , can the S trajectory be forced into transition to its equilibrium value. Figure 10 pictures dynamics of settling into equilibrium. Here, the natural capacity B∕ D = 10 7 , D = 1∕1500 , q B0 = 0.5 , B = 4 and other values are as in Sect. 6. If the average recovery period is 2.6 weeks, then 20 time units equal one year.

Compared to duration of an epidemic, dynamics of settling into equilibrium are slow, as is well known. (See, e.g., analysis of the linearized SIR model in Diekmann et al. 2013 .) For R B = 1.2 , as in Fig. 10 , the oscillation period is about 21 years and decay time constant about 140 years. For R B = 2 , they are about 9 and 90 years. The equilibrium I * is small compared to the epidemic I MAX : 78 people for R B = 1.2 and 235 for R B = 2 . Transitioning as in Fig. 10 , amplitude of oscillations in S and R depends on R B and imperceptibly (if at all, given numerical uncertainty) on size of the small incursion or its timing (within limits described earlier).

Computed values agree with Eq. (87) to within small fractions of a percent.

We review here main points of this article, summarize how the probabilistic model differs from other models in the literature and opine on advantages and disadvantages these differences confer.

Beginning with a small but nontrivial difference, the probabilistic model adopts a second-order gamma distribution for removal probability while most other models employ one of first order. The first-order distribution is attractive for its simplicity but is unrealistic in its time profile. In concert with the mass action law, the secondorder removal distribution predicts an expected maximum fraction of population infected at any one time, i MAX , significantly greater than does the first order, e.g., more than 30% for R 0 < 2.75 . This is consequential for predicting peak stress on medical facilities.

Our numerical studies indicate that the limit of the ratio of i MAX produced by the gamma distribution of order n to that of order one is 2n∕(n + 1) as R 0 goes to one.

Moving to a more significant difference, the probabilistic model incorporates two parameters versus one in both isolation and quarantine. We believe that modeling both strength ( q A0 and q B0 ) and speed ( A and B ) of these processes is framework for more realistic representation. These parameters are easy to interpret and potentially measurable, which should make the model useful in empirical research. When quarantine or combined intervention strength is modeled, we discover that intervention effectiveness, as measured by either s T∞ or i TMAX , is highly sensitive to small changes at a critical value of that strength. The critical value is about twothirds of the containment requirement derived from reproduction number. Does this accurately model the natural world? If so, knowledge of the phenomenon should be valuable.

Another major difference is that our model is founded on first principles of probability rather than mathematical simplicity. Coefficients such as 2 R B ∕( R + B ) that arise naturally in the probabilistic model have no counterpart in typical SIRlike models. Note that the coefficient combines parameters from different processes. The typical SIR-like model describes transitions with pairs of equations like ̇x i = − ij x i + ... and ̇x j = ij x i + ... , where the rate ij common to the pair is a parameter set by the analyst. Whether basis in probability is an advantage remains to be seen. We expect it to represent reality more accurately, but comparison with field data is required. A disadvantage: the ODE equations modeling quarantine are not time-invariant, which can be a hindrance, e.g., in stability analysis. We have modeled quarantine as time-driven, which we believe is more realistic than event-driven, e.g., removal only after tracing from contact with an infectious person.

An examination of the Reed-Frost theory of epidemics

On the spread of a disease with gamma distributed latent and infectious periods

An attempt at a new analysis of the mortality caused by smallpox and the advantages of inoculation to prevent it

An attempt at a new analysis of the mortality caused by smallpox and the advantages of inoculation to prevent it

Mathematical epidemiology: past, present and future

Mathematical models of isolation and quarantine

Modeling the worldwide spread of pandemic influenza: baseline case and containment interventions

When is quarantine a useful control strategy for emerging infectious diseases?

Some mathematical developments on the epidemic theory developed by Reed and Frost

Modeling the impact of quarantine during an outbreak of Ebola virus disease

Mathematical tools for understanding infectious disease dynamics

Endemic models with arbitrarily distributed periods of infection

Strategies for mitigating an influenza epidemic

Factors that make an infectious disease controllable

Mitigation strategies for pandemic influenza in the US

Modeling strategies for controlling SARS outbreaks

Epidemic disease in England-the evidence of viability and of persistence

Modeling target layered containment of an influenza epidemic in the US

How mathematical epidemiology became a field of biology

The mathematics of infectious disease

Effects of quarantine in six endemic models for infectious diseases

A contribution to the mathematical theory of epidemics

Analysis and numerical simulations of a stochastic SEIQR epidemic system with quarantine-adjusted incidence and imperfect vaccination. Special issue: Multiscale computational models for respiratory aerosol dynamics with medical applications

Mechanistic models of infectious disease and their impact on public health

Simulated pandemic influenza outbreaks in Chicago. Virginia Bioinformatics Institute

Generality of the final size formula for an epidemic of a newly invading infectious disease

Dynamics of an SEQIHRS epidemic model with media coverage, quarantine and isolation in a community with pre-existing immunity

Acknowledgement I thank reviewers for their constructive comments, and especially the reviewer who saved me from several mistakes in derivations and pointed me to an important reference.Funding Not applicable.

The author declares no conflict of interest.

We used reproduction numbers to determine requirements on q A0 , q B0 , A and B to contain epidemics of varying severity, confirmed them with ODE results and quantified dependence of intervention effectiveness on when interventions begin. We derived expressions for endemic disease equilibria when affected by isolation and found familiar formulas, with R B replacing R 0 . We examined dynamics of coming to an equilibrium and found that only special conditions permit it.In sum, assuming a homogeneous population we developed a model based on first principles of probability, with eight differential equations, six parameters and five reproduction numbers, and quantified important aspects of isolation and quarantine effectiveness in controlling an epidemic.