key: cord-0686737-hsorpcg1
authors: Vyasarayani, C.P.; Chatterjee, Anindya
title: New approximations, and policy implications, from a delayed dynamic model of a fast pandemic
date: 2020-08-25
journal: Physica D
DOI: 10.1016/j.physd.2020.132701
sha: 70014a9da05db32aef2349ec46342e7bf3b6b5f6
doc_id: 686737
cord_uid: hsorpcg1

We study an SEIQR (Susceptible-Exposed-Infectious-Quarantined-Recovered) model due to Young et al. (2019) for an infectious disease, with time delays for latency and an asymptomatic phase. For fast pandemics where nobody has prior immunity and everyone has immunity after recovery, the SEIQR model decouples into two nonlinear delay differential equations (DDEs) with five parameters. One parameter is set to unity by scaling time. The simple subcase of perfect quarantining and zero self-recovery before quarantine, with two free parameters, is examined first. The method of multiple scales yields a hyperbolic tangent solution; and a long-wave (short delay) approximation yields a first order ordinary differential equation (ODE). With imperfect quarantining and nonzero self-recovery, the long-wave approximation is a second order ODE. These three approximations each capture the full outbreak, from infinitesimal initiation to final saturation. Low-dimensional dynamics in the DDEs is demonstrated using a six state non-delayed reduced order model obtained by Galerkin projection. Numerical solutions from the reduced order model match the DDE over a range of parameter choices and initial conditions. Finally, stability analysis and numerics show how a well executed temporary phase of social distancing can reduce the total number of people affected. The reduction can be by as much as half for a weak pandemic, and is smaller but still substantial for stronger pandemics. An explicit formula for the greatest possible reduction is given.

Infection rate m Density of contacts β Net infectivity (βm) β n Baseline infectivity in a society β l Infectivity in the presence of social distancing or lockdown α

Rate of immunity loss γ Self-recovery rate p Probability of identifying and isolating an infected individual σ Asymptomatic and noninfectious time duration 

This work is partially motivated by the global pandemic of COVID-19. Understanding the dynamics of infectious diseases in a population can help in developing strategies to mitigate the spread [1, 2] . This paper presents new mathematical approximations and an asymptotic solution for a specific dynamic model for such infectious diseases. Some policy implications are discussed as well.

Mathematical models for the spread of disease have almost a century old history. In their seminal paper, Kermack and McKendrick [3] proposed a three-state model (popularly known as SIR) governing the evolution of susceptible (S), infected (I), and recovered (R) populations. In their model, the recovered population is assumed to have developed immunity against the infection. The model contains two free parameters, one for infection rate and one for recovery rate. The SIR model is widely used to predict the number of infected people in closed populations. The model has an analytical solution. Over time, the SIR model [4] has been modified to study infections where the recovered population can be reinfected (as with the common cold) and is known as the two-state SIS model. In the classic endemic model [5] , for diseases that are active over 10-20 years, information of new births and deaths are included. In another variant of the SIR model known as the four-state MSIR [6] model, passive immunity inherited by newborns from their mothers is included: for example, newborn babies can be immune to measles for 2 J o u r n a l P r e -p r o o f some time after their birth, but become susceptible later on. Other modifications have considered the effect of a carrier population [6] , which never recovers from the disease but is asymptomatic (relevant to, e.g., tuberculosis). Such people can again suffer from the disease later, or continue to infect others while remaining asymptomatic. In SEIR [7] , a four-state model, one of the states (E) represents the exposed population, infected but non-infectious. In the SEIQR model [8] , yet another state, representing a quarantined population, is added to the SEIR model. All the models discussed above, including SIR, SIS, MSIR, SEIR, and SEIQR are governed by nonlinear differential equations. More complicated partial differential equation models that include the effect of the age structure [4] of the population and vaccination history are also available [7] .

The models mentioned so far need not include time delays. However, the incubation, asymptomatic, and symptomatic phases of a disease can be incorporated as time delays in mathematical models. Including such delays in the differential equation models make them delay differential equations (DDEs), also known as retarded functional differential equations. Researchers have used DDEs [9, 10] to model the spread of infections like Zika [11] , HIV [12] , influenza [13] , and Hepatitis B [14] . Recently, a scalar DDE with one delay (time of recovery) based on a logistic model was used to study the spread of COVID-19 in Italy [15] . Some policy implications based on parameter studies and simulations were reported, but analytical progress on the delayed equation was limited.

We observe that the abovementioned models, with a few states, can be formally developed from underlying network models. Network epidemiological models with a large number of states [16, 17, 18, 19, 20, 21] follow fine details of the spatial and temporal spread of an infection. Their average or overall behavior can be described using what we might call lumped, compartmental, or continuum models.

In a recent article [22] , a network of compartmental models, with each model representing an age group was used to study the COVID-19 outbreak in India. Other works that have recently appeared, which study the progression of COVID-19 in different parts of the world, include [23, 24, 25, 26] . These works primarily deal with parameter identification, numerical simulations, and prediction of the number of cases with time, along with some policy implications.

Even with lumped models, if various effects [4] like incubation times, natural birth and death rates, prior immunity, and carrier states are included, then analytical solutions are usually unavailable. However, sometimes difficult problems can be solved approximately using asymptotic methods [27, 28] or related methods, and this paper presents such approximations for a slightly simplified case that is relevant to a fast-spreading pandemic.

In recent work that is directly related to our paper [9] , a five-state SEIQR system with delays has been developed as a continuum limit from a network model under quite general conditions. That system has then been examined as a generic model for infectious diseases. The model, even after simplification, has multiple parameters and coupled states, and so analytical progress is difficult. However, several interesting limits, steady states, and stability criteria have been presented [9] , along with supporting numerical results and policy implications.

Here we take up a simplified version of that model [9] . Our simplification is only that we ignore the possibility of some past sufferers of the disease eventually losing their immunity, and becoming vulnerable to infection all over again. With this simplification, significant new analytical approximations and insights are possible, and constitute the contribution of this paper. We note that Young et al. [9] say the following (their state of (1, 0, 0, 0, 0) is an uninfected population): complete evolution from infinitesimal infection to final saturation. These new approximations provide useful new analytical support and understanding that is not included within [9] . An overview of the paper is depicted in Fig. 1 . In section 2, the five-state SEIQR DDE model of Young et al. [9] is presented. Upon setting α = 0 (for a fast pandemic), E, Q, and R become slave variables. An SI DDE model with two-states is obtained.

A further simplified case of γ = 0 and p = 1 (zero self-recovery and certain quarantine) is considered in section 3. An asymptotic solution for weak growth using the method of multiple scales (MMS) has been obtained. For moderate growth, a long-wave approximation is considered, which converts the DDE into a first-order scalar ODE.

In section 4, the general case of γ > 0 and 0 < p ≤ 1 (nonzero self-recovery and probability of quarantine less than or equal to unity) is discussed. In this case, the long-wave approximation converts the DDE (SI model) into a second-order ODE. Excellent agreement is found between the solution of the ODE and that of the DDE for parameters associated with various diseases.

In section 5, the two-state SI model of section 2 is modified to account for time-varying β. A Galerkin approximation is used to convert the DDE model with time-varying β into six first-order ODEs. For both continuous and discontinuous variations in β, an excellent match is found between the solutions of the ODEs and the DDE. This match verifies the assumed low-order behavior of the DDE model, assumed throughout this paper.

In section 6, a policy of social distancing that involves piecewise constant variations in β in considered. In particular, a policy that uses a low value of β for an extended period, followed by a higher "normal" value, is shown to reduce the total fraction affected quite significantly (e.g., from 85% to 55% for some parameter choices). An explicit formula for the maximum reduction possible in this way is found. Notably, the reduction depends only on fraction affected, and not on any other system parameters. This provides a key policy implication and is the main practical contribution of the paper. As mentioned above, our mathematical model is essentially that of Young et al. [9] , with one of their small parameters set to zero. Figure 2 , adapted from their paper, shows the lumped model which is itself obtained by them from an underlying network model. In the model there are five subpopulations (actually population fractions) that add up to unity. The general governing equations [9] are:

In the above equations, the free parameters are interpreted as follows:β is the infection rate, m is the density of contacts, γ is the self-recovery rate, p is the probability of identifying and isolating an infected individual, and α is the rate of immunity loss. We have not introduced any simplifications of our own so far. We note, first, that E and Q in equations (2) and (4) are influenced by S and I along with their delayed values, but E and Q do not themselves influence S, I and R. In other words, E and Q are slave variables, and we henceforth ignore them.

In the three equations that remain, we can clearly absorb m intoβ or equivalently, writẽ

If social distancing is practiced, then we expect m to decrease and thus β to be lower, even thoughβ may remain the same. For this reason, in the later portions of this paper, we will consider time-varying β while holding other parameters constant. For a fast-spreading pandemic, we assume α = 0 for simplicity, which makes R a slave as well, and we need to only retain the equations for S and I. Finally, by choice of units of time, we can let σ = 1. This is equivalent to nondimensionalizing τ which has units of time, as well as γ and β which have units of 1/time. Our equations now arė

In the above, if β was constant rather than time-varying, then its t-dependence would be dropped. Note that, if β varies with time, its variation can be considered externally specified and not a part of the solution. Equations (7) and (8) make intuitive sense in a lumped-variable setting as follows. Equation (7) says that the instantaneous rate of new infections is proportional to how infectious the disease is (β), how much people are meeting each other (m(t)), how many uninfected people there are (S(t)) and how many infectious people are out in public (I(t)). Equation (8) says that the rate of change in the number of infectious people is equal to previously infected people just exiting the latency phase and entering the infectious phase, minus the rate at which people are recovering on their own, minus also the rate at which people displaying symptoms are being put into quarantine (these quarantined people are slightly diminished in number due to self-recovery, which is good; and due to some people not being quarantined, which is a system inefficiency).

We are interested in near-unity initial conditions for S(−∞) = 1 that lead to growth of the infection and eventual saturation. In particular, the net damage done by the disease is represented by 1 − S(∞). Infected individuals E(t) remain asymptomatic and non-infectious for a time duration σ. Subsequently, these individuals become infectious and enter population I(t), but remain asymptomatic for a time duration τ . Upon showing symptoms, they enter population Q(t) and are quarantined with probability p for a time κ, beyond which they infect nobody. Some infectious asymptomatic individuals may become non-infectious on their own, with rate γ. After quarantine, the cured population R(t) could in principle lose immunity at a small rate α, but we take α = 0 for a fast-spreading pandemic.

From now onward, until the start of section 5, we will consider β to be a constant parameter. For clarity, therefore, we write the governing equations out with constant β,

Let

where we are interested in the asymptotic initial condition P (−∞) = 0 as the limiting case of a tiny level of initial infection. Then equation (9) yields

which incorporates the initial condition of interest, namely S(−∞) = 1. Inserting S(t) from equation (12) into equation (10), we obtain

Integrating both sides with respect to time, and by defininḡ

where 1 −p is an integration constant chosen to match initial conditions at −∞. Thus, for constant β and with the approximation of α = 0 for a fast-spreading pandemic, equations (1) through (5) effectively reduce to the single nonlinear delay differential equation shown in (15) . It may be noted thatṖ = 0, i.e., P equals a constant, is allowed by equation (15) for P that satisfy

Equation (16) is satisfied by P = 0 for all parameter values. Additionally, for γ > 0, it has a single strictly positive root if

If γ = 0 (i.e., there is no self-recovery), and 0 ≤p < 1 (i.e., not everybody is quarantined), then for β > 0, P = 0 is the only equilibrium solution. This means if P increases from zero, it can grow without bound and S(∞) = 0, i.e., everybody in the population gets infected. The case of β = 0 is not interesting because the infection does not spread. Finally, if γ = 0 andp = 1, then equation (16) is identically satisfied for every constant P . Thus, we conclude that a simple yet interesting situation within equation (15) occurs when p = 1 (all infected people display symptoms and are quarantined) and γ = 0 (there is no recovery without displaying symptoms). We will first study this restricted case in some detail, because some analytical progress is possible that provides useful insights.

3 A simple subcase: p = 1 and γ = 0 For p = 1 and γ = 0, we have from equation (15),

As mentioned above, any constant P is an equilibrium, though possibly an unstable one.

In the initial stages P (t) is small, and equation (18) can be linearized tȯ

Equation (19) has infinitely many characteristic roots, and has oscillatory solutions which we must disallow because their decreasing portions require negative I(t). However, monotonic solutions exist as well, and we will examine them. The characteristic equation of equation (19) is

Among the infinitely many roots of the above equation, those with nonnegative real parts are of main interest because they lead to growth of P (t) from initial tiny values. We note that if the real part of λ is assumed nonnegative, then the magnitude of the right hand side is bounded by 2β. This means, for infinitesimal β, any right half plane roots of equation (20) are infinitesimal as well. We also note that λ = 0 is a root regardless of β.

For non-infinitesimal β, a criterion for real roots in the right half plane is easily found because λ = 0 is a root of equation (20) and the right hand side first increases and then decreases to zero as λ → ∞. These conditions imply that a nonzero positive root is assured if the slope of the right hand side at λ = 0 exceeds unity. That condition is βτ > 1.

Equation (21) suggests that the contagion may take hold if βτ > 1, and perhaps not if βτ < 1. Note that τ is fixed by biology but β can be lowered by practicing social distancing, which will be discussed towards the end of the paper. Incidentally, in Young et al. [9] , the instability condition including p ≤ 1 and γ > 0 is given as equation (17). For p = 1 and γ → 0, equation (17) easily reduces to equation (21).

For initial insight, we consider some numerical solutions of equation (18) obtained using Matlab's built-in solver dde23 with specified error tolerances of 10 −7 or better. The initial function used was

with a small and positive, and mentioned in the figure captions. Once P (t) is obtained (numerically or, later, analytically), we can calculate S(t) by using equation (12) . Results below will therefore be presented only in terms of the original variable S(t), since the impact of the disease is reflected by 

The initial function used for integrating equation (18) was P (t) = a 1 + t 1+τ , t ≤ 0. The two solutions in (b) are relatively time-shifted because the underlying system is autonomous, and one solution grows from smaller initial values; see main text for further discussion.

Two solutions for β = 1 and τ = 0.8 are shown in figure 3(a) . It is seen that not much decay occurs; and the solution in each case saturates at a value proportional to the magnitude of the initial function used. This is linear behavior (recall also the discussion following equation (21)).

In contrast, two solutions for β = 1 and τ = 1.2 are shown in figure 3 (b). Numerical solutions for two different initial functions show that the rapidly spreading phase of the disease and the saturation value of S are essentially identical, independent of the initial function. It is important to note that the relative time-shift between these two solutions is inconsequential. The underlying dynamical system is autonomous, and the use of asymptotic initial conditions at −∞ leaves a free parameter in the solution that allows time-translations (this will be explicitly clear in the analytical approximations later in the paper). The solution that starts from smaller initial values takes a little longer to climb to larger values, resulting in the time-shift.

Next, two solutions for β = 1 and β = 2, but with βτ = 1. (18) was

is the population fraction that remains unaffected, appears independent of β for βτ held fixed. Similar behavior is seen in the two curves in figure 4 (b) for βτ = 2 and identical initial conditions. Again, S(∞) appears independent of β for βτ held fixed.

Prompted by numerics, let

for some C to be determined. In the solution of interest, I(t) starts from infinitesimal values, and each individual who is a part of S(t) stays in the infectious state for exactly τ units of time. Since the total number of individuals thus affected is exactly 1 − S(∞), we can see that for the solution of interest 1

From equation (23), we obtain

Equation (25) needs to be solved numerically, but two limits are clear.

The other limit is for βτ 1, where C → βτ . Numerically, for βτ = 1.2, we find C = 0.376 or S(∞) = e −0.376 ≈ 0.687, which matches figure 3(a); and for βτ = 2, we find C = 1.594 or S(∞) = e −1.594 ≈ 0.203, which matches figure 3(b). The difference between even βτ = 1.2 and βτ = 2, though both are unstable, is large in terms of consequences for the population.

The above results indicate that if βτ > 1, then a solution that grows asymptotically from a zero value at minus infinity saturates at a value given by equation (26) . However, all P values are equilibrium values: it is therefore interesting to ask what the minimum value of P is for which the equilibrium is stable. We can find this by considering the corresponding limiting steady value of S to be S * , writing equation (10) for p = 1 and γ = 0 aṡ

and concluding that the required stability condition is (by adapting equation (21)) βτ S * < 1.

For βτ slightly exceeding unity, referring to equations (23) and (26), the upper limit of S * = 1 βτ implied by equation (27) corresponds to about half as many people being infected as would be if the asymptotic solution for constant β was allowed to run its course all the way from initiation to saturation. This situation will be clearly seen in the multiple scales solution for constant β and weak growth starting from zero at t = −∞. We now turn to that solution.

The foregoing results indicate that the P (t) solution starts asymptotically from zero at t = −∞ and grows monotonically if βτ > 1, but saturates at a small value when βτ slightly exceeds unity. We can develop an asymptotic solution for the case where

We will use the method of multiple scales [28, 29, 30, 27] . We rewrite equation (18) aṡ

and note that the = 0 case is on the stability boundary for P (t) = 0. We introduce three time scales for a second order expansion,

10 J o u r n a l P r e -p r o o f Journal Pre-proof think of P (t) as P (T 0 , T 1 , T 2 ) with a slight abuse of notation, and observe that the time derivative is to be interpreted asṖ

Further, a delayed quantity such as P (t − ∆) is to be interpreted as

where due to the smallness of , Taylor series expansions in can be used for the second and third arguments, but not the first argument, i.e.,

Finally, P itself is to be expanded as

where higher order terms in the expansion would require retention of still slower time scales. This much is routine, and yields an equation of the form (using the symbolic computation software Maple; note that the leading order term is O( )):

where two long expressions have been written simply as L 2 and L 3 (details omitted for brevity). At O( ) we have

for which we adopt the solution (based on previous observations; and also upon rejecting fast-varying exponential decaying terms [29, 30] )

i.e., P 0 is a constant on the fast or T 0 time scale. Inserting equation (37) into equation (35), we obtain at O( 2 ),

with terms containing A(T 1 , T 2 ) canceling each other out at this order. This means A(T 1 , T 2 ) remains indeterminate at this order, and we are free to choose

because it adds nothing new to the already retained P 0 . Inserting the above into equation (35) , we obtain at O( 3 ),

J o u r n a l P r e -p r o o f where L 3 is a long expression independent of T 0 and containing the function A(T 1 , T 2 ) along with its T 1 -derivatives only. Now, appealing to the required boundedness of P 2 (which corresponds to removal of secular terms, and can also be viewed as a choice that allows the approximation to stay valid for a longer time), we insist that L 3 = 0. Further, since T 2 derivatives of A do not appear, we set A back to a function of T 1 alone (no contradiction up to this order). In this way, we finally obtain

where we note that τ is a positive parameter, A is a function of T 1 , and primes denote T 1 -derivatives. The above differential equation is to be solved as a function of T 1 , with the initial condition A(−∞) = 0. The solution turns out to be

where c 0 is an indeterminate constant that allows time-shifting (recall figure 3 (b) and the related discussion). Inserting T 1 = t, and then (recall equation (28)

we finally obtain the leading order approximation for the entire solution, for βτ slightly greater than unity, as

where c 1 is an undetermined constant. The above solution saturates, as t → ∞, at

which matches equation (26) upon noting that we can replace C β with Cτ with no errors introduced at leading order. At the same order of approximation, we can also write

whence

where we have used the '∼' notation because this is an asymptotic approximation. A numerical example is given in figure 8 for βτ = 1.025 (β = 1 and τ = 1.025) and βτ = 1.05 (β = 1 and τ = 1.05). The match is good, but deteriorates for larger values of βτ . It may be noted that, for such weak growth where only a small fraction of the total population gets infected before the pandemic runs its course, by equation (46),

In comparison (recall equation (27)), S could in principle be stable at a value as high as

i.e., the uncontrolled and constant-β solution infects about twice as many people as seems strictly necessary; we will return to this in section 5. We will next develop a long-wave approximation that performs better for somewhat larger βτ . (18) is P (t) = a 1 + t 1+τ , t ≤ 0. The arbitrary time-shift c 1 of the multiple scales solution (equation (47)) is chosen here to obtain a good visual match.

The advantage of using an asymptotic method like the method of multiple scales is that we have formal validity as → 0. We are reassured that the solution we are seeking does in fact exist, and has approximately the shape obtained as the leading order approximation. However, in the present case, for somewhat larger the solution is not very accurate; moreover, proceeding to higher orders leads to long expressions that seem difficult to simplify usefully. Therefore, encouraged by our asymptotic solution, we now develop a more informal but more accurate approximation. In particular, we try a long-wave (LW) approximation as follows.

We suppose that there is a "long" scale (technically a time scale, for this problem) which we shall call L, such that

Let

Now equation (18) 

Expanding the above in a series for large L, retaining terms up to O(L −2 ), and solving forP , we obtain

We note that on the right hand side the second term is dominant because it contains the large parameter L, while the first term does not. The left hand side has the highest derivative and cannot be dropped 13 J o u r n a l P r e -p r o o f without changing the order of the differential equation. For these reasons, a further approximation to equation (53) isP

Equation (54) is integrable, and yieldŝ

where C 0 is an integration constant. The initial condition of interest isP = 0 andP = 0 as ξ → −∞, whence we obtainP

.

We now return to equation (53), where we can insert equation (56) into the non-dominant term as an approximation. A simple way to do it is to replace the first term on the right hand side with an approximation that is integrable, i.e.,

This givesP

The above is integrable, and after enforcing the initial conditionP = 0 andP = 0 as ξ → −∞, setting L = 1 (it was a bookkeeping parameter all along, helping us keep track of which effects are big and which are small), and dropping the hat, we obtain the informal long-wave approximation

We mention that equation (59) is indeed a significant simplification, because it replaces a DDE (infinite-dimensional phase space) with a first order ODE (one-dimensional phase space). It also applies to a specific solution, all the way from infinitesimal initiation to final saturation. The freedom in a single initial condition that it offers is equivalent to simply the arbitrary time-shift already noted in the multiple scales solution. While it cannot be solved explicitly in closed form, its solution can be formally expressed in implicit form using an indefinite integral (omitted for brevity).

The qualitative behavior of the approximation in equation (59) may be checked easily for small positive P by linearizing the right hand side to obtaiṅ

which is consistent with the earlier result that growth requires βτ > 1 (inequality (21) ). In figure 6 , numerical solutions for βτ somewhat greater than unity are shown. We consider βτ (18) is P (t) = a 1 + t 1+τ , t ≤ 0. The initial condition P (0) is used to integrate the long-wave equation (59). The coefficient c 1 is used in the multiple scales solution (see equation (47)).

Finally, a direct analytical comparison with the multiple-scales solution can be made by expanding the right hand side of equation (59) in a power series for small P , retaining upto quadratic terms. That equation can be solved in closed form, and gives a hyperbolic tangent solution as well, which initially looks a little different from the multiple scales solution:

However, noting that the multiple scales solution was for βτ close to unity, if we replace (2 − βτ ) with 2 − 1 = 1, and replace βτ (τ + 2) with 1 · 1 β + 2 , then we recover equation (46). This is not surprising because for βτ slightly greater than unity, the long-wave approximation is asymptotic as well. The approximation is only informal (as opposed to asymptotic) when we use it for arbitrary values of β and τ , with βτ not close to unity.

This concludes our study of the p = 1 and λ = 0 subcase, which is interesting both because it is a reasonable limit and because it permits any constant P as an equilibrium. The latter is not true for general parameter values, to which we turn next. Encouraged by the simplicity of the long-wave approximation as opposed to the multiple scales solution, we try only the former.

Proceeding with the same ansatz as equation (50), expanding up to second order, and setting L to unity, we obtainP

where

A few things may be noted here. First, assuming that τ andp are not too small, we assume µ > 0 in equation (63). Second, we were able to use a trick for the p = 1 and γ = 0 case to reduce the order of the approximating long-wave differential equation; for those parameter values, the last term on the right hand side in equation (62) becomes zero. Here, we have had to retain the second order ordinary differential equation. Third, if there is a nonzero positive P for which equation (62) has an equilibrium solution, then that P must satisfyp

which is equivalent to equation (16) . This means the steady state P in this general case will be exactly correct, unlike the previous subcase of p = 1 and γ = 0. However, note that this unique limiting P is for the specific solution that starts asymptotically at zero, at t = −∞. Since we now have a second order differential equation in the long-wave approximation, trying to approximate what is purportedly a single solution, we have a minor dilemma in interpreting the two degrees of freedom available in initial conditions for equation (62). One of those corresponds to an arbitrary time-shift, as noted earlier. That leaves one other initial condition. Since our derivation has not identified what this initial condition should be, we expect that it should be inconsequential. This does turn out to be the case, as seen next. For small P and smallṖ , equation (62) can be linearized to givë

where C 1 and C 2 are given by

For stability of P = 0 we must have C 1 > 0 and C 2 > 0, but this solution is of little interest because the infection does not grow. If C 1 < 0 and C 2 > 0, then growing oscillations are predicted. Such solutions are non-physical for our application, because P must be monotonic. For monotonic growth with a single unstable growth direction, we require C 1 > 0 and C 2 < 0. We assume that µ > 0 (recall equation (62) and equation 63) accordingly our conditions become

Equation (68) matches the loss of stability condition of Young et al. [9] as well as our own equation (17) . In other words, the condition for a solution growing from zero exactly matches the condition for the 16 J o u r n a l P r e -p r o o f existence of another equilibrium for strictly positive P . This condition, though obtained here from the long-wave approximation, matches the theoretical result exactly because it is near this stability boundary that the long-wave approximation is asymptotic. Note that if C 1 > 0 and C 2 < 0 in equation (65), then the two characteristic roots are real: one is positive and one is negative. The positive root leads to the growing solution, which is of primary interest in this paper. The negative root absorbs the apparently free initial condition, contributes an exponentially decaying term that dies soon, and has no influence on the growing solution provided the initial conditions are sufficiently early in the outbreak. A numerical example is shown in figure 7 for the case of C 1 > 0 and C 2 < 0. The parameter p = 0.98 implies that the probability of detecting infected individuals is not perfect. Further, γ = 0.1 models a small fraction of the infectious population recovering without being quarantined. The Long-wave solution almost perfectly matches the solution from numerical integration of the DDE (equation (18)). The phase portrait shown in figure 7 (b) indicates that the stability properties of the solution of interest (shown in red) are correctly inherited by the long-wave approximation. Nearby solutions are attracted towards the same asymptotic solution starting from infinitesimal initial infection and growing monotonically towards final saturation.

For this specific kind of solution, the long-wave approximation gives a good match to the DDE solution as long as µ > 0, C 1 > 0 and C 2 < 0.

In figure 8 , we compare the solution of the long-wave ODE (equation (62)) and that of the DDE (equation (15)) for parameters corresponding to six different diseases. These parameters are reported in Table 1 and are taken from references [9, 31, 32, 33, 34, 35, 36] . In figure 8 , the long-wave solution captures the disease progression accurately in all cases. We now consider more general dynamics, starting from less restricted initial conditions, and allowing for a time-varying infection rate β.

As a part of its public health policy, based on observed spread of the disease, a government may prescribe temporarily greater social distancing. Social distancing effectively lowers β. We can then use equations (7) and (8), rewritten here for readability (incorporatingp from equation (14)):

It is worthwhile to verify that the system has low dimensional behavior. To this end, in this section we present a six-state Galerkin approximation for equations (70) and (71) using Legendre polynomials as basis functions. The ODEs from the Galerkin approximation are shown below (see the appendix for a derivation and [37, 38] for details). Writing ν = 1 + τ , α 2 = 1 − 2 ν and α 3 = 3 2 1 − 2 ν 2 − 1 2 , the reduced order model iṡ Figure 8 : Comparison between numerical solution (of the DDE) and long-wave solutions for parameter sets corresponding to six different diseases (a) H1N1 (Brazil) [9] (b) Ebola [31] (c) Spanish Flu 1917 [33] (d) SARS [33] (e) Hepatitis A [33] , and (f) COVID-19 [34, 35] . Initial conditions and other parameter values used to generate the results are shown in each figure saparately.

J o u r n a l P r e -p r o o f

In the above approximation S(t) = η 1 (t) + η 2 (t) + η 3 (t) and I(t) = η 4 (t) + η 5 (t) + η 6 (t). Here, β(t) is a known function of time and the state variables do not have delays. Note that we approximate the dynamical system here, and not a specific solution as we did with the long wave approximation. The accuracy of the reduced order model is demonstrated in figure 9 . For each of several sets of parameters, we find an excellent match between numerical solutions of the DDEs and the Galerkin based ODEs. Both continuous and discontinuous β's are considered. The reduced order model shows that our DDEs (equations (70) and (71)), though formally infinite-dimensional systems, are effectively finite-dimensional. The remaining dynamics consists of rapidly decaying components that are soon inconsequential.

In the next section, we present some policy implications of time-varying β.

Since social distancing lowers β, we consider the implications of time-varying β in the simple case where β is set to a low value β l early in the pandemic, held at that value for a sufficiently long time, and then set back to its higher, normal value β n . We emphasize that β n is a matter of culture and lifestyle, which people may not like to change permanently; and use of β l < β n is a temporary measure. We will show that the implications can be significant in terms of reducing the total number of people affected. Moreover, a simple analytical formula for the obtainable benefits will be found. Consider a country with a constant β = β n under normal living conditions. If the pandemic runs its course without any social distancing measures, the solutions of equations (70) and (71) approach S(t) = S n and I(t) = 0. To study the stability of these solutions, we substitute S(t) = S n + ξ and leave I unaltered, and drop the subscript n from β, to obtain from equations (70) and (71)

The stability of equation (79) determines the stability of the system. Substituting I = e λt into equation (79), and setting λ = 0 at the bifurcation point, the condition for stability is seen to be

From equation (64), the steady-state value S n satisfies (note the n subscript on β)

It is seen above that if γ > 0 butp = 1, then S n = 1. This means that, for large p and small τ , S n will not be much smaller than unity. However, in practice, γτ is not very small, and p is not very large. As an example representative of COVID-19 (γ = 0.08, p = 0.2, τ = 10, and β n = 0.1962, andp = pe −γτ ), we find that S n = 0.15 from equation (81). The possible advantages of social distancing are now seen from the inequality 80. Substituting β n = 0.1962, S n = 0.15, along with the abovementioned values forp and γ, we find the inequality 80 is met by a relatively large margin: 0.0268 < 0.08.

In fact, for the same β = β n , stability could potentially be achieved by a significantly larger steady value of S(t), say 

J o u r n a l P r e -p r o o f which works out to 0.45 for the parameter values being considered in this example. In other words, the number of unaffected people could potentially be three times larger, and the number of affected people could be less by about one third (0.55 as opposed to 0.85). The policy implication is that if we hold β at a lower value β l and let the pandemic run its course and stabilize at S(t) = S L , then subsequently raising β back to β n will yield a stable state, and further infections will not grow. Stated another way, if β is held at a constant value for all time, the pandemic progresses beyond the point of first herd immunity to a state where more people are affected than necessary. However, with a temporary phase of social distancing, the controlled pandemic can in principle grow to precisely the level where herd immunity just becomes effective for the normal level of interactions, at which point social distancing can be lifted with no further penalties. In this way, while preserving our long term levels of social interactions, we can achieve the smallest possible numbers of people affected.

The precise value of β l needed is easy to compute for the present model. We first use our normal parameters in equation (82) to compute the best achievable outcome, S L . The corresponding β = β l is found using equation (64), which here becomes

which for the parameter values of the present example yields β l = 0.1279. The simple social distancing policy recommendation now is

where social distancing is imposed at time T o , early in the pandemic when very few people have been affected; and T c is understood to be sufficiently late such that an intermediate equilibrium at S L is established. Subsequently, β is set back to β n and life proceeds as usual.

In figure 10 , we show the evolution of the pandemic for three cases: (i) β(t) = β n = 0.1962 (constant for all time), (ii) β(t) = β l = 0.1279 (held constant for all time), and (iii) β(t) given by equation (84) with T 0 = 50 and T c = 550. In case (i), 85% of the population is affected. In case (ii), only 55% of the population is affected. In case (iii), too, only 55% is affected even though β is later increased to the level of case (i).

The effectiveness of this variable-β policy can be expressed using the ratio ρ < 1, defined by

If ρ is significantly less than unity, then the benefits of adopting the above optimal social distancing policy are high. A simple expression for ρ can be obtained as follows. First, equation (82) is solved for β n , which is inserted into equation (81). Simple manipulations then yield

The above formula is in terms of S n , which is formally a quantity that depends on model parameters β, τ , γ and p. However, in terms of interpretation, the dependence on S n alone is both interesting and useful, because it shows the answer is dependent only on the severity of the pandemic. In particular, for very small S n , ρ ≈ 0.5. However, even for S n = 0.15 (i.e., 85% of the population would normally be affected), we find that ρ = 0.65, or about two thirds. In the popular press, the idea is commonly expressed that social distancing spreads the infected people out by flattening the curve, and this gives the medical services of a country more time to respond. Our result shows that the total number affected can be significantly reduced as well. 

In this paper we have taken up a recently presented SEIQR model with delays. For a fast-spreading pandemic, loss of immunity of previously infected and cured people may reasonably be ignored. Under that simplification, the SEIQR model decouples so that only the S and I population equations need to be tackled. It is known for this model that, for fixed parameter values in the unstable regime, an outbreak can occur. An initially small infected population can grow, and a significant portion of the original population can be affected.

We have first studied this model in some detail, seeking useful approximate solutions. For a weakly growing outbreak that affects a small proportion of the total population, and under a further simplification that neglects self-recovery and assumes perfect quarantining, the method of multiple scales yields an analytical expression for the complete progression of the outbreak, from infinitesimal initiation to final saturation. For moderate growth rates, a long wave approximation for the same parameters provides a nonlinear first order ODE for the same progression. With imperfect quarantining and nonzero self-recovery the long wave approximation for the full progression of the outbreak is given by a second order ODE. Finally, although the underlying DDE system is technically infinite dimensional, we have shown that a six-state Galerkin-based reduced order model for the system does an excellent job of capturing a wide range of solutions, i.e., the dynamics is effectively low-dimensional.

Subsequently, we have examined the implications of policy-induced social distancing, incorporated in our model as a time-varying infection rate β(t). Interestingly and promisingly, we have found that an extended period of social distancing, imposed early in the outbreak, followed by an eventual relaxation to usual levels of interaction, can significantly lower the total numbers infected without losing stability of the final state. In the limit of weak growth, the number of infected people is cut in half. For faster growth, the reduction is a smaller but still significant. We have obtained a simple analytical formula for the reduction possible.

The above policy implications seem simple and robust. The intuitive key to understanding this reduction caused by social distancing lies in stability under fresh, but small, infection. Here, stability implies that with a small infected population, the outbreak will not grow very much (recall figure 3(a) versus 3(b)). Under identical conditions, a larger infected population could cause the outbreak to grow: 23 J o u r n a l P r e -p r o o f the assumption is that once the infected numbers are contained, a fresh large influx of infected people will be avoided. If β is set to β l with social distancing, the outbreak saturates at a relatively high S L . Subsequently, assuming no large influx of infected people, β can be increased back to the normal β n , and the outbreak does not grow significantly further.

An introduction to infectious disease modelling

Mathematical structures of epidemic systems

A contribution to the mathematical theory of epidemics

Mathematical models in population biology and epidemiology

Three basic epidemiological models

The basic epidemiology models: models, expressions for R 0 , parameter estimation, and applications

The mathematics of infectious diseases

An SEIQR model for childhood diseases

Consequences of delays and imperfect implementation of isolation in epidemic control

Numerical modelling in biosciences using delay differential equations

A Fractional-Order Model for Zika Virus Infection with Multiple Delays

Mathematical analysis of delay differential equation models of HIV-1 infection

Dynamics of a delay differential equation model of Hepatitis B virus infection

Solvable delay model for epidemic spreading: the case of Covid-19 in Italy

Epidemic model with isolation in multilayer networks

Six Susceptible-Infected-Susceptible models on scale-free networks

Efficiency of prompt quarantine measures on a susceptible-infectedremoved model in networks

Rapid decay in the relative efficiency of quarantine to halt epidemics in networks

Epigrass: a tool to study disease spread in complex networks

Networks and epidemic models

Age-structured impact of social distancing on the COVID-19 epidemic in India

Data analysis and modeling of the evolution of COVID-19 in Brazil

2020 A Mathematical Description of the Dynamics of Coronavirus Disease (COVID-19): A Case Study of Brazil

Mechanistic-statistical SIR modelling for early estimation of the actual number of cases and mortality rate from COVID-19

2020 Management strategies in a SEIR model of COVID 19 community spread

Perturbation methods

Multiple Scales and Singular Perturbation Methods

Multiple scales without center manifold reductions for delay differential equations near Hopf bifurcations

Order reduction of retarded nonlinear systems-the method of multiple scales versus center-manifold reduction

Estimating the Reproduction Number of Ebola Virus (EBOV) During the 2014 Outbreak in West Africa

Comparing nonpharmaceutical interventions for containing emerging epidemics

The reproductive number of COVID-19 is higher compared to SARS coronavirus

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Galerkin Projections for Delay Differential Equations

Galerkin approximations for stability of delay differential equations with time periodic coefficients

Here we outline our Galerkin projection calculation. Readers interested in the theoretical background may see, e.g., the so-called tau method of imposing boundary conditions in [39] .The initial functions for equations (70) and (71) are assumed to be S(t) = U 1 (t), −ν ≤ t ≤ 0 and I(t) = U 2 (t), −ν ≤ t ≤ 0. Define y 1 (s, t) = S(t+s) and y 2 (s, t) = I(t+s). Equations (70) and (71) along with their history functions can be equivalently posed as the following partial differential equations with time dependent boundary conditionsNow, we assume a solution for y 1 (s, t) and y 2 (s, t) as followsThe basis functions φ 1 (s) = 1, φ 2 (s) = 1 + 2s ν , and φ 3 (s) = 3 2 1 + 2s The inner products with φ 3 (s) are not taken. Instead, we substitute (90) and (91) in the boundary conditions (88) and (89). There, we have y 1 (0, t) = η 1 + η 2 + η 3 , y 2 (0, t) = η 4 + η 5 + η 6 , y 1 (−1, t) = η 1 + φ 2 (−1)η 2 + φ 3 (−1)η 3 , y 1 (−ν, t) = η 1 − η 2 + η 3 , y 2 (−1, t) = η 4 + φ 2 (−1)η 5 + φ 3 (−1)η 6 , and y 2 (−ν, t) = η 4 − η 5 + η 6 . Equations (88) and (89) becomėgiving us six ODEs for the six states, equivalent to equations (72)-(77). The initial conditions for our ODEs can be obtained from history functions as η k (0) = 2k−1 ν 0 −ν U 1 (s)φ k (s)ds, k = 1, 2, 3 and η r (0) = 2(r−3)−1 ν 0 −ν U 2 (s)φ r−3 (s)ds, r = 4, 5, 6.

J o u r n a l P r e -p r o o f