New approximations, and policy implications, from a delayed dynamic model of a fast pandemic

1 Introduction

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 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):

“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.

2 Mathematical Model

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:

• Download figure
• Open in new tab

Figure 1:

The SEIQR model with delays. Healthy individuals S(t) are infected with rate constant β. 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 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, write

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 are

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 defining we obtain where 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 , 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 (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 and , 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.

3.2 Initial Numerical Observations

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⁻⁷ 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)).

• Download figure
• Open in new tab

Figure 2:

Two solutions each for (a) βτ = 0.8 < 1 and (b) βτ = 1.2 > 1. The initial function used for integrating equation (18) was . 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.

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 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.2 in each case, and with identical initial functions (a = 0.0001) are shown in figure 3(a). The limiting or saturation value S(∞), which 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(∞) appears independent of β for βτ held fixed.

• Download figure
• Open in new tab

Figure 3:

Two solutions for (a) βτ = 1.2 > 1 and (b) βτ = 2 > 1. The parameters β and τ are individually varied. The initial function used for integrating equation (18) was .

3.5 Multiple Scales Solution for Weak Growth

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) as 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 (t) as P (T₀, T₁, T₂) 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 𝒪(ϵ)): where two long expressions have been written simply as L₂ and L₃ (details omitted for brevity). At 𝒪(ϵ) 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₀ is a constant on the fast or T₀ time scale. Inserting equation (37) into equation (35), we obtain at 𝒪(ϵ²), with terms containing A(T₁, T₂) canceling each other out at this order. This means A(T₁, T₂) remains indeterminate at this order, and we are free to choose because it adds nothing new to the already retained P₀. Inserting the above into equation (35), we obtain at 𝒪(ϵ³), where L₃ is a long expression independent of T₀ and containing the function A(T₁, T₂) along with its T₁-derivatives only. Now, appealing to the required boundedness of P₂ (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₃ = 0. Further, since T₂ derivatives of A do not appear, we set A back to a function of T₁ 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₁, and primes denote T₁-derivatives. The above differential equation is to be solved as a function of T₁, with the initial condition A(−∞) = 0. The solution turns out to be where c₀ is an indeterminate constant that allows time-shifting (recall figure 2(b) and the related discussion). Inserting T₁ = ϵ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₁ is an undetermined constant. The above solution saturates, as t → ∞, at which matches equation (26) upon noting that we can 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 βτ.

• Download figure
• Open in new tab

Figure 4:

Comparison between numerical solution and asymptotic (method of multiple scales, or MMS) solution for (a) β = 1, τ = 1.025, a = 1 × 10⁻⁵ (numerical) and c₁ = 310 (MMS); (b) β = 1, τ = 1.05, a = 1 × 10⁻⁵ (numerical) and c₁ = 200.5 (MMS). The initial condition used for integrating equation (18) is . The arbitrary time-shift c₁ of the multiple scales solution (equation (47)) is chosen here to obtain a good visual match.

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 βτ.

3.6 Long-Wave Approximation for Moderate Growth

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) becomes

Expanding the above in a series for large L, retaining terms up to O(L⁻²), and solving for , 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) is

Equation (54) is integrable, and yields where C₀ is an integration constant. The initial condition of interest is and as ξ → −∞, whence we obtain

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 gives

The above is integrable, and after enforcing the initial condition and 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 obtain 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.

• Download figure
• Open in new tab

Figure 5:

Comparison between numerical, multiple scales (MMS), and long-wave solutions for (a) β = 1, τ = 1.15, a = 1 × 10⁻⁵ (numerical), c₁ = 93 (MMS), and P (0) = 0.87 × 10⁻⁴ (long-wave); (b) β = 1, τ = 1.25, a = 1 × 10⁻⁵ (numerical), c₁ = 66 (MMS), an d P (0) = 0.78 × 10⁻⁴ (long-wave). The initial condition used for integrating equation (18) is . The initial condition P (0) is used to integrate the long-wave equation (59). The coefficient c₁ 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 solutio n was for βτ close to unity, if we replace (2 − βτ) with 2 − 1 = 1, and replace βτ (τ + 2) with , 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.

4 Long-Wave Approximation for General Parameter Values

We now take up equation (15), reproduced below:

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

Where

A few things may be noted here. First, assuming that τ and 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 satisfy 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 give

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).

• Download figure
• Open in new tab

Figure 6:

Comparison between numerical and long-wave solutions for (a) β = 1, τ = 1.15, a = 1 × 10⁻⁵ (numerical), and P (0) = 16.9 × 10⁻⁷ (long-wave); (b) β = 1, τ = 1.25, a = 1 × 10⁻⁵ (numerical), and P (0) = 4 × 10⁻⁷ (long-wave). The initial condition used for integrating equation (18) is . The initial condition P (0) reported in the figure and are used to integrate the long-wave equation (62).

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 β.

5 Time-Varying β and Related Policy Implications

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 (incorporating 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 and the reduced order model is

In the above approximation S(t) = η₁(t) + η₂(t) + η₃(t) and I(t) = η₄(t) + η₅(t) + η₆(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.

• Download figure
• Open in new tab

Figure 7:

Comparison between numerical and Galerkin solutions for various parameters. H(t − c) = 1 if t > c, and is zero otherwise. Initia l functio ns used for equations (67) and (68) are and , t ≤ 0. Initial conditions for equations (69)-(74) are fitted by Galerkin projection (see the appendix).

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.

• Download figure
• Open in new tab

Figure 8:

Role of time-varying β. For β(t) = 1.04 held fixed, 7.6% of the population gets infected. With β(t) = 1.04 − 0.02H(t − 300) + 0.02H(t − 2500), 3.9% of the population gets infected and the outbreak proceeds no further. This is practically the same as having β = 1.02 for all time. When the switch occurs later, with β(t) = 1.04 − 0.02H(t− 500)+0.02H(t − 2500), a slightly larger fraction of the population gets infected. Initial conditions used for equations (67) and (68) were and .

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.

• Download figure
• Open in new tab

Figure 9:

Influence of T_c in β(t) = 1.04 − 0.02H(t − T_c) + 0.02H(t − T_o) on the percentage of the total infected population at final stable steady state. The parameters used are the same as in figure 8, except for T_c and T_o.

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.

6 Concluding Discussion

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.

Appendix

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₁(t), −ν ≤ t ≤ 0 and I(t) = U₂(t), −ν ≤ t ≤ 0. Define y₁(s, t) = S(t+s) and y₂(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 conditions

Now, we assume a solution for y₁(s, t) and y₂(s, t) as follows

The basis functions and are shifted Legendre polynomials defined on the domain −ν ≤ s ≤ 0. Substitute (79) into (75) and (80) into (76). Premultiplying each equation with ϕ₁(s) and then by ϕ₂(s) and integrating over the domain −ν ≤ s ≤ 0 each time, we obtain

The inner products with ϕ₃(s) are not taken. Instead, we substitute (79) and (80) in the boundary conditions (77) and (78). There, we have y₁(0, t) = η₁ + η₂ + η₃, y₂(0, t) = η₄ + η₅ + η₆, y₁(−1, t) = η₁ + ϕ₂(−1)η₂ + ϕ₃(−1)η₃, y₁(−ν, t) = η₁ − η₂ + η₃, y₂(−1, t) = η₄ + ϕ₂(−1)η₅ + ϕ₃(−1)η₆, and y₂(−ν, t) = η₄ − η₅ + η₆. Equations (77) and (78) become 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 and .