key: cord-0228554-8rv6i0hp authors: Vyasarayani, C. P.; Chatterjee, Anindya title: New approximations, and policy implications, from a delayed dynamic model of a fast pandemic date: 2020-04-08 journal: nan DOI: nan sha: 4a2974ce8fe71bcb034e757c05867e9ba71fbaf5 doc_id: 228554 cord_uid: 8rv6i0hp We study an SEIQR (Susceptible-Exposed-Infectious-Quarantined-Recovered) model 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 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 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 correctly executed time-varying social distancing, within the present model, can cut the number of affected people by almost half. Alternatively, faster detection followed by near-certain quarantining can potentially be even more effective. 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 "The fact that orbits starting from near (1, 0, 0, 0, 0) will approach an endemic equilibrium is difficult to prove; this part of our prediction is supported by numerical simulations." In this context, we will present a new asymptotic multiple-scales solution for weak growth in a special case, and two informal long-wave approximations for moderate growth, all three solutions describing the complete evolution from infinitesimal infection to final saturation. These new approximations provide useful new analytical support and understanding that is not included within [9] . We will also discuss transient solutions and the role of initial conditions, in the context of a time-varying infection rate, using a six-state Galerkin approximation based reduced-order model obtained from the original delayed equations. Since social distancing affects the infection rate, we will demonstrate how, within the present model, properly executed social distancing can cut the total number of affected people by almost one half. A stronger case for early detection and quarantine will emerge simultaneously. As mentioned above, our mathematical model is essentially that of Young et al. [9] , with one of their small parameters set to zero. Figure 1 , 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ė 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. 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(∞). 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ḡ we 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. 5 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 1 − S(∞). Two solutions for β = 1 and τ = 0.8 are shown in figure 2 (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 2 (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 (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 3 (b) for βτ = 2 and identical initial conditions. Again, S(∞) 7 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 2(a); and for βτ = 2, we find C = 1.594 or S(∞) = e −1.594 ≈ 0.203, which matches figure 2(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)) For βτ slightly exceeding unity, referring to equations (23) and (26), the upper limit of 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. 1 An analogy may help explain this trick. Imagine a room where a finite but large number of people, say M people, enter over a long period of time and at a variable rate. The number of people in the room, N (t), is an arbitrary nonnegative function of time. Each person stays in the room for exactly τ units of time, and then leaves. Clearly, 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, 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 ), 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 2(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 4 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 Now equation (18) becomesP 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 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 (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)). 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 5 , numerical solutions for βτ somewhat greater than unity are shown. We consider βτ = 1.15 (see figure 5(a)) and βτ = 1.25 (see figure 5 (b)). The numerically obtained long-wave solutions match well with full numerical solutions of the original DDE; in fact, they match significantly better than the multiple scales solution. 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 , 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. We now take up equation (15) , reproduced below: Proceeding with the same ansatz as equation (50) 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ë For stability, coefficients of bothṖ and P above need to be negative. If the coefficient ofṖ is positive while the coefficient of P remains negative, then growing oscillations are predicted; such solutions are non-physical for our application, because P must be monotonic. For monotonic growth, it is the coefficient of P that must be positive. This leads to the condition for existence of such solutions, which matches the (loss of) stability condition of Young et al. [9] as well as our equation (17). In other words, the condition for a solution growing from zero is exactly the same condition for the existence of another equilibrium for strictly positive P ; and 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. Finally, when the coefficient of P in equation (65) is indeed positive, then the two characteristic roots are real: one is positive and one is negative. The positive root leads to the growing solution, as above. 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 6 for two cases: βτ = 1.15 and βτ = 1.25. 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. In the figure, the long-wave solution perfectly matches the final saturation state obtained from numerical integration of the DDE (equation (18)), due to the equivalence of equations (64) and (16) . We have so far concentrated on approximations for a specific solution: one that starts from infinitesimal infection levels, changes monotonically but slowly over a long time and accelerates over a relatively short time, to finally saturate at a finite value as the time goes to infinity. We now consider more general dynamics, starting from less restricted initial conditions, and allowing for a time-varying infection rate β. As public health policy based on observed spread of the disease, a government may prescribe temporarily greater social distancing, effectively lowering β. We can then use equations (7) and (8), rewritten here for readability (incorporatingp from equation (14)): We first present a six-state Galerkin approximation for equations (67) and (68) using Legendre polynomials as basis functions. The ODEs from the Galerkin approximation are shown below (see the appendix for a brief derivation and refer to [31, 32] for details). Writing ν = 1 + τ , α 2 = 1 − 2 ν and α 3 = 3 2 1 − 2 ν 2 − 1 2 , the reduced order model iṡ In the above approximation S(t) = η 1 (t) + η 2 (t) + η 3 (t) and I(t) = η 4 (t) + η 5 (t) + η 6 (t). Here, β(t) has delays but the state variables do not. 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 shown in figure 7 . 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 (67) and (68)), though formally infinite-dimensional systems, are effectively finite-dimensional. The remaining dynamics consists of rapidly decaying components that are soon inconsequential. Having demonstrated that the dynamics is effectively low-dimensional, we can have greater faith in simple numerical studies that suggest policy implications. We now turn to such a policy implication. The issue was anticipated in section 3, and we now discuss it in detail. We suppose that normal living entails some specific β, and the public cannot indefinitely maintain high social distancing, i.e., significantly lower β. Yet it may be possible to lower β early in the outbreak, and then go back to the normal β later, when it is safe. The benefits are illustrated using simulations in figure 8 , where the chosen τ , γ and p correspond to a critical β = 1 (recall equation (66)). Now, suppose that normally β = 1.04. The disease will spread, and saturate at S(∞) ≈ 0.92, as per equation (48). Yet, for the same β = 1.04, a larger uninfected population of S * could be stable (equation (49)). That S * , in principle, could be reached using an artificially low β ≈ 1.02. If we change β from 1.04 to 1.02 early in the outbreak, and hold β = 1.02 until a steady state is reached, then finally returning to β = 1.04 could yield a stable solution. The outbreak would be arrested, and the total number of affected people would be cut almost in half. In figure 8 , β is switched from 1.04 to 1.02 at two different instants T c , for two different simulations, and then switched back to 1.04 later. In each such switched case, although finally β = 1.04, the steady value S(∞) corresponds to β = 1.02. In contrast, if β = 1.04 had been held throughout, almost twice as many people would have been infected. A further numerical study is presented in figure 9 . Here we vary both T c as well as T o , the time of switching back to β = 1.04. The parameters of figure 8 are used except for T c and T o . We see that the percentage of people affected can be cut almost in half for sufficiently low T c provided T o is large enough. If T o is made smaller, then there is a special value of T c when S(∞) is highest for the chosen T o , but the gains may be suboptimal. For higher β, more people must get infected (immune) before stability is achieved, and the benefit obtained is proportionately a little lower. For example, for the same τ , γ and p, if β = 1.1 is held fixed, then the final uninfected population is 82.4%. In contrast, upon switching down to β = 1.05 for an extended period before returning to β = 1.1, the final uninfected population is 90.2%, i.e., the number of infected people decreases by 46%. The tradeoffs between T c and T o are similar to those observed in figure 9 , except that the dynamics is about two times faster (graphical results omitted for reasons of space). The reader may anticipate that similar benefits can be obtained by lowering the time to quarantine, τ , guided by the instability criterion of inequality (21) . However, with p < 1 and γ > 0, the actual instability criterion is inequality (66). Here, for example, if p = 0.9 and γ is tiny, then the stability condition is almost independent of τ . However, if p is close to unity and γ is not too small (as in the parameters used in the foregoing examples), then the benefits of reducing τ can be significant. To demonstrate with a numerical example, the steady value S(∞) can be computed from equation (64) using With β = 1, p = 0.98 and γ = 0.1, we obtain table 1. It is seen that relatively small reductions in τ achieve significant reductions in the net population affected by the disease. However, if p = 0.92 (i.e., the chance of escaping quarantine is four times greater) and γ remains the same, then instability occurs at τ = 0.22 (about four times smaller), and S(∞) is less sensitive to τ as well (numerical examples omitted for brevity). We conclude that in circumstances where p can be kept high (i.e., if public institutions are strong and detection followed by quarantine is nearly certain), and where the self-recovery rate γ is not extremely small, the system can benefit significantly from even modest reductions in the detection time τ . In such situations, research directed toward earlier detection may yield substantial benefits, and it may even be unnecessary to engage in extreme social distancing. On the other hand, if there is under-reporting and imperfect quarantining, then impractically large reductions in τ may be needed to compensate, and strict social distancing may be more effective. 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 little smaller. Additionally, if the probability of an infected person being detected and quarantined is high, and the self-recovery rate not too small, then perhaps even stronger benefits can be obtained by slightly reducing the time until quarantine, τ . The above policy implications seem simple and robust. We emphasize that the benefit is not merely in lowering the rate at which people get infected, but also in the total number of people infected by the end of the outbreak. 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 2(a) versus 2(b)). Under identical conditions, a larger infected population could cause the outbreak to grow: the assumption is that once the infected numbers are contained, a fresh large influx of infected people will be avoided. If β is lowered with social distancing, the outbreak saturates at a high S(∞), and the infected population goes to near-zero values. Subsequently, under the assumption of no subsequent large influx of infected people, β can be increased within the stability boundary, and the outbreak does not grow significantly further. In contrast, the benefits from reducing the time to quarantine, τ , require greater sustained institutional alertness but may be stronger. In closing, we must acknowledge that in any lumped model of the sort we study here, spatial variations in parameters and infected population densities are not modeled. Such lumped models are averaged models. Thus, it is not really clear at a detailed spatial level what it means to reduce the average infection rate β by, say, 2%. If the average person engages in social distancing, benefits will be seen on average, although there could still be localized outbreaks within pockets where people engage in riskier behavior. In constrast, however, the speed and probability of being detected and quarantined is less up to individual members of the public, and more in the hands of public institutions. Such institutions are amenable to tighter quality measures. So, from the viewpoint of reliability, we believe that large scale testing and near-certain isolation or quarantining can be critically useful in containing pandemics like COVID-19. 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 A delay differential model for pandemic influenza with antiviral treatment 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-infected-removed 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 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 [33] .The initial functions for equations (67) and (68) 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 (67) and (68) 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 (79) and (80) in the boundary conditions (77) and (78). 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 (77) and (78) becomėgiving us six ODEs for the six states, equivalent to equations (69)-(74). 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.