Modelling a pandemic with asymptomatic patients, impact of lockdown and herd immunity, with applications to SARS-CoV-2.

The SARS-CoV-2 is a type of coronavirus that has caused the pandemic known as the Coronavirus Disease of 2019, or COVID-19. In traditional epidemiological models such as SEIR (Susceptible, Exposed, Infected, Removed), the exposed group E does not infect the susceptible group S. A distinguishing feature of COVID-19 is that, unlike with previous viral diseases, there is a distinct "asymptomatic" group A, which does not show any symptoms, but can nevertheless infect others, at the same rate as infected symptomatic patients. This situation is captured in a model known as SAIR (Susceptible, Asymptomatic, Infected, Removed), introduced in Robinson and Stillianakis (2013). The dynamical behavior of the SAIR model is quite different from that of the SEIR model. In this paper, we use Lyapunov theory to establish the global asymptotic stabililty of the SAIR model, both without and with vital dynamics. Then we develop compartmental SAIR models to cater to the migration of population across geographic regions, and once again establish global asymptotic stability. Next, we go beyond long-term asymptotic analysis and present methods for estimating the parameters in the SAIR model. We apply these estimation methods to data from several countries including India, and demonstrate that the predicted trajectories of the disease closely match actual data. We show that "herd immunity" (defined as the time when the number of infected persons is maximum) can be achieved when the total of infected, symptomatic and asymptomatic persons is as low as 25% of the population. Previous estimates are typically 50% or higher. We also conclude that "lockdown" as a way of greatly reducing inter-personal contact has been very effective in checking the progress of the disease.

Introduction

Background

The mathematical modelling of the spread of epidemics has a long history, stretching back over several centuries. The “modern” approach to the modelling of epidemics can be said to have begun with Kermack and McKendrick (1927), which first enunciated the principle that infected persons pass on the disease to susceptible persons at a rate proportional to the number of contacts between the two groups. Over the years, various refinements of the basic model have been proposed. The literature on disease modelling is truly enormous. Indeed, a survey paper (Hethcote, 2000) published in 2000 already had more than 200 references. Today it would be many times that number. Therefore it would be futile to attempt a summary of the entire topic of epidemic modelling. Rather, in Section 2 we limit ourselves to highlighting those aspects of epidemic modelling that are broadly common to most existing models, and why these models are not adequate to study the latest health-related challenge, namely the onset of the COVID-19 pandemic.

Traditionally, epidemiologcal models have grouped people into two, three or four groups, usually denoted by Susceptible (S), Exposed (E), Infected (I), and Removed (R). Note that many authors use the symbol R to denote “recovered.” Note that in the epidemiology literature, the phrase “compartment” is widely used instead of “groups.” Moreover, we prefer to use “groups” because, later in the paper, we study the impact of migration on the spread of a pandemic using compartmental models. However, in this paper we wish to make a distinction between those who recover and are immune to reinfection, and those who die from the disease. This distinction becomes particularly important when we introduce births and deaths due to natural causes into the model. In traditional models, contact between a member of the infected group I and another person belonging to the susceptible group S leads to the latter person becoming infected with a certain probability. Depending on the model, the susceptible person either becomes infected straightaway (the SIR model), or enters an intermediate stage called Exposed (E) (SEIR model). In the latter scenario, it is assumed that contact between persons belonging to the E and S groups does not lead to fresh infections, because members of the E group do not carry a sufficient viral load to infect others through contact.

However, one of the characteristic features of the coronavirus pandemic is that many of the persons who contract the disease are “asymptomatic,” or belong to the group A. A recent paper (Oran & Topol, 2020) collates several publicly available data points and states that “Asymptomatic persons seem to account for approximately 40% to 45% of SARS-CoV-2 infections, and that they can transmit the virus to others for an extended period of time, perhaps longer than 14 days.” Other references estimate the fraction of asymptomatic patients to be more than 50% at times (Mizumoto, Kagaya, Zarebski, & Chowell, 2020), and as high as 75% (Day, 2020). For this reason, asymptomatic patients remain “hidden” and cannot be identified except through testing the entire population, which is clearly impractical.

Moreover, asymptomatic patients (A) differ from exposed patients (E) in one important respect. Unlike in traditional epidemiological models, contact between a person in the A group and another in the S group does lead to the latter getting infected, with a certain probability. In addition, as in other models, contact between a person in the I group and another in the S group also leads to the latter getting infected, with a similar probability. To the best of the authors’ knowledge, the first paper to formulate and analyze a model that captures this phenomenon is Robinson and Stilianakis (2013). However, the analysis of the model in Robinson and Stilianakis (2013) is not so complete as is currently available for the SEIR model. As shown in Korobeinikov and Maini (2004); Korobeinikov and Wake (2002) and reviewed in subsequent sections of this paper, the global stability properties of the SEIR model are well understood, where those of the SAIR model are still being studied. The authors of Robinson and Stilianakis (2013) did not give a name to their model. In the present paper, we adopt the model of Robinson and Stilianakis (2013) and refer to it as the SAIR (Susceptible, Asymptomatic, Infected, Removed) model. Then we carry out a complete analysis of the behavior of this model, both with and without vital dynamics, on a par with what is currently known about the SEIR model. As shown in later sections, the two dynamical models are quite different, as are the solutions.

Organization and contributions of the paper

The paper is organized as follows: We begin by reviewing two classical models, namely the SIR and the SEIR models, and analyze the stability of these models using Lyapunov stability theory. This analysis closely follows Korobeinikov and Maini (2004); Korobeinikov and Wake (2002). Then we carry out a complete analysis of the SAIR (Susceptible, Asymptomatic, Infected, Removed) model, both without and with vital dynamics. In order to establish the stability properties of the SAIR model, we extend the classical Krasovskii-LaSalle theory of Lyapunov stability, from the case where the Lyapunov function V is positive definite, to the case where V is only positive semidefinite. This extension is of independent interest. Then we introduce the idea of “compartmental” models, wherein there are multiple compartments within a society, each of them having its own SAIR groups, with different levels of interaction and other epidemiological parameters. Note that in the epidemiology literature, the phrase “compartment” is often used to denote the various groups S, E, I, R. However, in the dynamical systems literature, “compartmental models” refer to collections of individual dynamical systems that interact with each other. We too use the phrase “compartment” in this sense, and refer to S, A, I, R as “groups.” We are able to extend the earlier stability analysis of SAIR model to a two-compartment model without vital dynamics. However, extending the SAIR model to a multi-compartment model with vital dynamics remains an open problem.

In the remainder of the paper, we first present methods for estimating the parameters in the SAIR model based on the evolution of the pandemic. Then we present the outcomes of applying our theories to actual data from the COVID-19 pandemic in eight countries from around the world, including India. Then we focus on the progress of the pandemic in Delhi, one of the “hotspots” in India. We model the effect of a “lockdown” whereby contacts between persons is severely limited, and show that quantitative predictions based on our models faithfully reproduce actually observed data. Finally, we examine whether the notion of “herd immunity” which has been propagated by some persons is real or not. We conclude that herd immunity is not only real, but is also achieved at far lower levels of community infection than was thought earlier. We conclude the paper by discussing several interesting problems that merit the attention of the research community.

Review of the SIR and SEIR models

As mentioned in the Introduction, the literature on epidemiological modelling is vast, and there is no point is even attempting a comprehensive review. Rather, we review two classical models known as the SIR and SEIR models respectively, and analyze their stability properties using Lyapunov theory. The SIR model forms the point of departure for the SAIR model, which provides a more realistic model for COVID-19, compared to the SEIR model.

The SIR model

SIR Model without vital dynamics

In the SIR model, the population is divided into three groups, denoted as S (Susceptible), I (Infected), and R (Removed). Note that many authors use R to denote “Recovered.” However, in our model, the group R also includes those who die from the disease. Also, it is assumed that the total population size is constant, so that S, I, R represent the fraction of the population within each group. ThereforeIn the simplest models, births and deaths (due to natural causes) during the course of the epidemic are not taken into account. In more detailed models, births and deaths (due to natural causes) are included, and these are called “vital dynamics” in Hethcote (2000). However, it is assumed that the births and deaths balance exactly, so that the overall population size remains constant. The introduction of vital dynamics significantly changes the dynamical behavior of the model. Specifically, without vital dynamics, the SIR model exhibits a continuum of equilibria, whereas with vital dynamics, the SIR model has one unstable and one attractive equilibrium (under suitable conditions). The assumption of constant population size can be removed as in Korobeinikov and Wake (2002) by replacing the quantity S by another quantity that those authors call P, which is less intuitive. It turns out that there is no difference between the behavior of the dynamical system whether the population size remains constant or not. Therefore we limit ourselves to the case of constant population.

In the absence of vital dynamics, the equations that govern the SIR model are Fig. 1 contains a flowchart of these equations. As expected, we have that . Therefore we can ignore anyone of the three equations and focus only on the other two. Most authors ignore R and studywhere β, γ > 0 are parameters of the disease under study. The logic behind (2) is as follows: When a person in group S makes contact with a person in group I, the former gets infected with a likelihood of β. Left to themselves, the persons in group I move to the group R (get removed) at a rate of γ. The ratio β/γ is referred to as the basic reproduction ratio, and is denoted by σ.1 As we shall see below, if σ ≤ 1, the pandemic does not take off, and dies down steadily, whereas if σ > 1, the pandemic initially grows before dying down.

Fig. 1

Flowchart of the SIR model without vital dynamics.

To make these ideas precise, let us analyze the dynamics in (2). Because we are ignoring R, this dynamical system evolves over the simplexIt can be seen that any point (S, 0) where S ∈ [0, 1] is an equilibrium of the system (2). Therefore there is a continuum of equilibria. This is consistent with the observation that the Jacobian matrix of the right side of (2) around any such equilibrium is a singular matrix.2 Therefore analysis methods based on linearlization around an equilibrium do not apply to this system.

A thorough analysis of this equation is carried out in Hethcote (1976); see Equation (2.5) and thereafter. Because we will be making use of these ideas in our SAIR model, we briefly reproduce the relevant details.

(See ( , Theorem 2.2).) Consider the system (2) starting at an initial condition (I 0, S 0). If σS 0 ≤ 1, then I(t)↓0 as t → ∞. If σS 0 > 1, then I(t) increases at first and then decreases to 0, while S(t)↓S ∞, where S ∞ is the unique solution in (0, 1/σ) of

By dividing the second equation in (2) by the first, we getIf we make the reasonable assumption that so that then the solution of (4) iswhere and . The behavior of the solutions is completely captured by the constant σ. When t → ∞, it is evident that and which in turn implies from (2) that . Substituting this into (5) readily gives (3). □

Next we discuss the concept of “herd immunity.” Though this term is introduced in Topley and Wilson (1923) (i.e., even earlier than the SIR model), the term did not have a precise definition nor analysis until the publication of Dietz (1975); Smith (1970), with the latter paper being more mathematical. A good summary of the evolution of the concept is found in Fine, Eames, and Heymann (2011). In diverse publications, the term “herd immunity” has been used to mean two apparently different things, namely: (i) the value of S(t) at which the level of infection I(t) is maximum, and (ii) . It is not clear whether the research community realized that, in the SIR model, both definitions are equivalent. Hence we formally state and prove this.

For a given σ, S 0 with σS 0 > 1, and the maximum value of I(t) occurs when and is given by where At this time instant, we have that

Note that I(t) assumes its maximum value when . Therefore, if I(t) ≠ 0, then if and only if . Substituting this value of S into (5) giveswhich is (6) after collecting terms. Now the fact that for all t implies that, when we have (7). □

From Theorem 2, it is clear that the number of infections is maximum precisely when ; thus, as indicated above, both usages of “herd immunity” are consistent. There is yet another usage of the phrase “herd immunity” that is found in the literature, namely: the level of (that is, the “immune” or non-susceptible population) at which the infection level I(t) begins to decrease monotonically to zero. For the sake of clarity, let us denote this value by H (to suggest “Herd”). The analysis in Theorem 2 shows thatThis particular formula for the herd immunity level is quite widely used in epidemiology. For instance, in Angulo, Finelli, and Swerdlow (2020), the authors suggest that, in the context of COVID-19, the total of recovered and currently infected persons must exceed 55% before various containment strategies (such as lockdown) can be relaxed. This number is arrived at by taking for COVID-19 and applying (8).

In any infection, it is reasonable to assume that I 0 ≪ 1. With this assumption, one can rewrite (6) asIf I 0 ≪ 1, thenSubstituting this estimate into (6) gives (9). From (9), it can be seen that as the basic reproduction ratio σ increases, the constant C(σ) increases towards 1, while the slope of the straight line approximation decreases. Fig. 2 shows the dependence of on I 0 for .

Fig. 2

Dependence of on I0 for various σ.

We conclude our discussion of herd immunity by revisiting an old formula (Dietz, 1975, Eq. (23)), which states (in current notation) thatThe logic behind the above equation is that it is difficult if not impossible to infer the basic reproduction ratio σ while an epidemic is in progress. However, once it has died down so that S ∞ is known, the above formula can be used to estimate σ. Now let us relate this formula to the current analysis. It is easy to rearrange (3) asIf I 0 ≪ 1 so that S 0 ≈ 1, this expression is nearly the same as (10).

The above analysis shows clearly that, in the absence of vital dynamics (natural births and deaths), every (S, 0) with S ∈ [0, 1] is an equilibrium. Moreover, for every (S 0, I 0) such that the corresponding solution converges to (S ∞, 0) (and of course ), where S ∞ is given by (3). This point can be quite far away from (S 0, I 0). In this sense, no equilibrium is stable in the sense of Lyapunov. However, the entire set of equilibria is globally attractive, in the sense that as t → ∞, the distance from (S(t), I(t)) to the set approaches zero. Note that, actually, S(t) → S ∞ which is in (0, 1/σ). This conclusion is arrived at by directly solving the SIR equations.

SIR Model with vital dynamics

Now we introduce vital dynamics, that is, births and deaths due to natural causes. The equations in this case are as follows:In this equation, α is the rate of birth as well as death. It is assumed that all newborns enter the S group, and in each of the groups S, I, R, people die at the rate of α. Note that now consists of two components: A term γI which corresponds to the infected persons being removed from the pool of the infected, either through immunity or death, and another term αR which corresponds to death due to natural causes. Note that births and deaths balance, so that we still haveSo once again we are able to ignore one of the three equations, namely for R. Fig. 3 presents a flowchart of the SIR model with vital dynamics.

Fig. 3

Flowchart of the SIR model with vital dynamics.

Some papers such as Korobeinikov and Wake (2002) examine the case where births and deaths do not balance. In this case, the original variables (S, I, R) can no longer be viewed as fractions, because the underlying population is itself changing with time. In Korobeinikov and Wake (2002), a transformation is presented whereby S is replaced by another variable P so that (P, I, R) satisfy . But P is not the same as S. For these reasons, we stay with the assumption of balanced births and deaths. After we complete the stability analysis for this situation, the reader may compare the results with those in Korobeinikov and Wake (2002) and see that unbalanced births and deaths do not change the qualitative behavior of the dynamical system.

Now let us determine the equilibria of the system (11). Unlike with (2) which has a continuum of equilibria, in this case there are only one or two isolated equilibria. Thus the introduction of vital dynamics actually simplifies the dynamics, as we shall see. One equilibrium, which is referred to as the nonendemic equilibrium, is . Another one (S*, I*) is known as the endemic equilibrium, and corresponds to I* ≠ 0, that is, there is a persistent level of infection at the equilibrium. Setting and dividing by I* (which is assumed to be nonzero) givesNote that this computation is meaningful only if S* ≤ 1, or equivalently . For simplicity we ignore the case and assume in the sequel that . In this case the quantityis called the “effective” reproduction rate and is more than one. Also, the introduction of vital dynamics decreases σ, because . Next, setting and substituting for S* givesNow let us return to the stability of the system (11). The analysis of this system has evolved over the years. Initially, the stability of these equilibria is analyzed using linearization in Hethcote (1976). It is shown that the nonendemic equilibrium is unstable, while the endemic equilibrium is (locally) asymptotically stable. In Li and Muldowney (1995), it is shown that the nonendemic equilibrium is in fact globally attractive, except for the other equilibrium. This is shown by establishing that the system does not have any limit cycles, by applying the Poincaré-Bendixson theorem (Vidyasagar, 2002, Theorem 3.3.22). However, that approach is inherently limited to two-dimensional systems. In Mena-Lorca and Hethcote (1992), a Lyapunov function is proposed, but it is not very elegant nor easy to analyze. The best solution to date is the Lyapunov function proposed in Korobeinikov and Wake (2002), which we present and study next. But before that we recall a well-known result often paraphrased as “the arithmetic mean is no smaller than the geometric mean.”

Suppose are positive numbers. Then with equality if and only if all x with equality if and only if for all i.

The proof is easy and is therefore omitted.

Now we establish the attractivity of the endemic equilibrium.

Define the set Suppose that . Then all trajectories starting in converge to the endemic equilibrium (S*, I*) as t → ∞.

We follow Korobeinikov and Wake (2002) and propose the Lyapunov function candidatedefined over the region . Now Therefore the Hessian of V is positive definite and V is strictly convex over . Moreover, its only stationary point is which is perforce the unique global minimum.3 Next, we computeThe computation is tedious but routine, and can be simplified by writingWhen the dust clears, we getNow apply Lemma 1 with and . Then it follows from (16) thatTherefore the function is positive definite, and everywhere. Moreover vanishes on the set {(S*, I): I ∈ [0, 1]}. Now a routine application of the Krasovskii-LaSalle invariance theorem (see e.g. Vidyasagar, 2002, Theorem 5.3.77) shows that the endemic equilibrium (S*, I*) is globally attractive over the set . □

The SEIR model

The SEIR model differs from the SIR model in that there is an additional group, known as Exposed (E). These are people whose viral load is not sufficient to infect anyone through contact.

SEIR Model without vital dynamics

The SEIR model without vital dynamics is described byThe above equations mean that when a person from group S comes into contact with a person from group I, then the former becomes “exposed” at a rate of β. Note that the transition is out of group S but to group E, and not to group I. The persons in group E become infected at a rate γ, and move to group I. Finally, people in group I move to group R at a rate of δ. Note that the transition of people is strictly sequential in the order S → E → I → R. A provision to move directly from group E to group R could be added with more burdensome notation. Fig. 4 contains a flowchart of these equations.

Fig. 4

Flowchart of the SEIR model without vital dynamics.

Note that there is no term of the form ES in the above equations. Therefore, contact between a susceptible person and an exposed person does not have any consequences. This is precisely the difference between previous diseases to which the SEIR model has been applied, and COVID-19. As before, we can ignore the equation for and focus on the other three. However, the equation for is useful to infer that at any equilibrium, we must have . It is easy to see that the set of equilibria of the system consists of all vectors of the form (S, 0, 0), S ∈ [0, 1].

Now let us introduce vital dynamics into the system. This model below is analogous to (11). As before, the model consists of adding a birth term α to and subtracting a multiple by α in all terms. This givesWe can streamline the equations by defining new constants which turns the above equations intoNote that a slightly more general model is used in Korobeinikov and Maini (2004). A flowchart of the SEIR model with vital dynamics is shown in Figure 5

Fig. 5

Flowchart of the SEIR model with vital dynamics.

As is the case in the SIR model with vital dynamics, there are now just two isolated equilibria, one with and one with I* ≠ 0, which are called the nonendemic and the endemic equilibria, respectively. The nonendemic equilibrium is . To compute the endemic equilibrium, we proceed as follows: Suppose I* ≠ 0. ThenThe expression for I* can be rearranged asNote that an endemic equilibrium exists only if θϕ ≤ βγ. In this case one can define the basic reproduction ratio as (βγ)/(θϕ). We will revisit this issue again towards the end of the next section. As before we ignore the possibility that and assume that θϕ < βγ.

Define the set Then, whenever the trajectory (S(t), E(t), I(t)) of the system (19) converges to the endemic equilibrium (S*, E*, I*).

In analogy with (17), we propose the Lyapunov function candidateAs before, this function is strictly convex over and has a global minimum at (S*, E*, I*). After some character-building computations, it can be shown thatNote that the second term is always nonpositive, by Lemma 1. As for the first term, note that the product of the three fractions is one. Hence, by Lemma 1, unless all three fractions are equal to one, that is, unless this term is negative. So is negative definite, and all solutions starting in converge to the endemic equilibrium. □

Note that, in contrast to the SIR model with vital dynamics, this time is negative definite, and there is no need to invoke the Krasovskii-LaSalle theorem.

Before concluding this section, we point out that the SIR and SEIR models discussed here contain what is usually referred to as “bilinear” incidence. Specifically, the infection resulting from contact between members of S and I is of the form βIS. It is possible to replace this by more general “nonlinear” incidence models, which is done in Liu, Hethcote, and Levin (1987). In that paper, the term βIS is replaced by βI, and some results are proved. Comprehensive results are proved in Korobeinikov and Maini (2004). The authors show that if p ≤ 1, then the endemic equilibrium has the region (in the case of SIR) and (in the case of SEIR) in its domain of attraction, by constructing a Lyapunov function that generalizes (21), and showing that is negative definite. It is noteworthy that the value of q > 0 is immaterial to the conclusions. A still more general class of incidences is studied in Korobeinikov and Maini (2005), and a complete stability analysis is carried out. Another related approach is found in O’Regan, Kelly, Korobeinikov, O’Callaghan, and Pokrovskii (2010), which studies the SIR model. In this paper, contrary to the usual practice, the equation for is ignored and the equations for are used to construct a Lyapunov function.

The SAIR model

In the preceding section we reviewed a couple of standard epidemiological models, and established their stability. In the SEIR model, it is assumed that contact between the Susceptible (S) and the Exposed (E) groups does not lead to fresh infections. However, one of the distinguishing features of COVID-19 is the presence of a large fraction of “asymptomatic” people who can, and do, infect susceptible people. To capture this phenomenon, we divide the population into four groups, denoted as Susceptible (S), Asymptomatic (A), Infected (I), and Removed (R). We examine the SAIR model both with and without vital dynamics, and analyze their stability properties. Note that the SAIR model is introduced in Robinson and Stilianakis (2013). However, to date it does not seem to have attracted much attention. Perhaps this is because, until the advent of COVID-19, there haven’t been any major epidemics with large asymptomatic patient populations. Some of the papers that cite Robinson and Stilianakis (2013) are Chisholm, Campbell, and Wu (2018); Mathur and Narayan (2018); Saad-Roya, Wingreena, Levinc, and Grenfell (2020). These papers do not go beyond the original analysis in Robinson and Stilianakis (2013). However, (Chisholm et al., 2018) has a good analysis on why ignoring asymptomatic patients can lead the policy-makers to under-estimate the severity of the pandemic.

We begin by stating and proving a result of independent interest.

An extension of the Krasovskii-Lasalle invariance theory

In this section we present an extension of the classical Krasovskii-LaSalle invariance theory, which is useful in establishing the stability properties of various epidemiological models studied here. The extension covers the case where a Lyapunov function candidate is only positive semidefinite, and not positive definite as is commonly assumed. The traditional Krasovskii-LaSalle theory is found in several sources, e.g., (Vidyasagar, 2002, Section 5.3), Actually the extension is implicitly contained in existing proofs of the theorems given in Vidyasagar (2002, Section 5.3), but was apparently not noticed earlier – perhaps because the need for it did not arise. However, in the analysis of the SAIR model, this extension is very useful.

The set-up is the standard one. We consider the differential equationwhere for some integer n, and the vector field f is continuously differentiable (so that the equation has a unique solution at least locally). Suppose satisfiesso that x* is an equilibrium. The objective is to derive conditions under which x* is “globally” attractive, in the sense that all solutions starting in some suitably large set (not just an immediate neighborhood of x*) converge to x*.

For this purpose we identify a Lyapunov function candidate that is continuously differentiable, and satisfies the following conditions:

There is a constant c such that the level setis compact (closed and bounded).

V satisfies

satisfies

Assumption [A2] is the key differentiator, because the usual assumption isWith these preliminaries, we can now state the extension:

Define and let Ω denote the largest invariant set of the system (22) contained in . Then x(t) → Ω as t → ∞, whenever . In particular, if {x*} is the only invariant set of (22) contained in then x(t) → x* as t → ∞, whenever .

The proof is fairly straightforward and follows well-established lines. Since for all it follows that the set is invariant for the flow of the differential Eq. (22). Since this level set is compact, the limit point set of each trajectory is nonempty, and the trajectory converges to its limit set; see Vidyasagar (2002, Definition 5.2.27) and Vidyasagar (2002, Lemma 5.2.34). Next, we refer to Vidyasagar (2002, Lemma 5.3.71). Though this lemma is stated for the case where V is positive definite, a perusal of the proof shows that it holds even if V satisfies only (24). Hence the desired conclusion follows. □

A related result can be found in Khalil (2002, Theorem 4.4, Corollary 4.1).

An SAIR model without vital dynamics

As mentioned earlier, the distinguishing feature of the COVID-19 pandemic is the presence of a large number of asymptomatic patients, who do not manifest any external symptoms, but are still capable of infecting susceptible persons. To capture this phenomenon, the following model is presented in Robinson and Stilianakis (2013):4

In the above model, S, A, I, and R denote the susceptible, asymptomatic, infected, and removed populations respectively. Interactions between A and S lead to the person from S moving to A at the rate of β, while interactions between I and S lead to the person from S moving to A at the rate of β. Note that persons from S move only to A and do not move directly to group I. The persons in group A move to the group R at the rate γ, and to the group I at the rate δ. Finally, persons in group I move to group R at the rate γ. Fig. 6 contains a flowchart of the above equations.

Fig. 6

Flowchart of the SAIR model without vital dynamics.

An alternative to the above model is to permit some fraction of the term to enter the group I directly, instead of passing through A as an intermediate stage. There does not appear to be any biological justification for this. Another possibility, which is totally unrealistic, is to combine into one group, and assume that A and I each make up a fixed fraction of the total. This would be just the SIR model “in disguise” with playing the role of I. From (26) it can be seen that initially the growth would be in group A which leads to growth in group I later on. This temporal behavior seems to tally with actually observed evolution of the pandemic. So we use the model in (26) throughout. However, we point out a variant of the SEIR model that could be interpreted as an SAIR model. In this SEIR-variant, (18) is modified towhere ϵ is a “small” number denoting secondary infections due to interactions between E and S; see van den Driessche and Watmough (2008, Section 6.4.2). As soon as the model includes the possibility that interactions between E and S lead to infections, this is the SAIR model, but for the restriction that ϵ is “small.” Having said that, we reiterate that Robinson and Stilianakis (2013) is apparently the first paper to use the acronym SAIR, and to formulate the model as in (26) without insisting that β.

It is easy to verify that the setis an invariant set of (26), and that the set of equilibria is {(S, 0, 0), S ∈ [0, 1]}.

Define For the system (26) , we have that

To analyze the stability of this system, we introduce the Lyapunov function candidateIt might be mentioned that the above function does not look very “traditional.” Nevertheless, it is positive definite over and has its global minimum at (0,0,0). NowHence on . Moreover, the set where vanishes is precisely and is an invariant set of the system (26). (In fact it is just the set of equilibria.) Now the desired conclusion follows from Theorem 5. □

Note that the above approach can also be applied to the SIR model using the Lyapunov function to show that I(t) → 0 as t → ∞. However, the simple nature of the SIR model allows us to draw much more detailed conclusions as in Hethcote (1976), and presented here as Theorem 1.

Now let us impose the simplifying assumptionsin (26). This leads toWe refer to this model as the simplified SAIR model, to distinguish it from the more general SAIR model of (26). Fig. 7 contains a flowchart of the simplified SAIR model.

Fig. 7

Flowchart of the simplified SAIR model without vital dynamics.

There are several noteworthy points about the simplified SAIR model, and these are discussed before proceeding to an analysis of this model.

It is assumed that the likelihood of fresh infection is the same, whether the contact is between A and S, or between I and S. Note that in Robinson and Stilianakis (2013), it is not assumed that these two rates are the same. This paper considerably predates the emergence of COVID-19. After the onset of the COVID-19 pandemic, several papers in the literature study “viral shedding” by both asymptomatic and infected patients, and conclude that there is no discernible difference between the two; see for example He, Lau, Wu, and Deng (2020); Li, Pei, and Chen (2020); Liu, Yan, and Wan (2020); Wölfel, Corman, and Guggemos (2020). Therefore, in the simplified SAIR model, we assume that the rate of infection due to S and A interactions is the same as that due to S and I interactions, and the same constant β is used to multiply the terms AS and IS.

It is assumed that, irrespective of the cause of infection, all infected persons enter only group A; this is similar to the assumption in the SEIR model that all infected persons enter group E. Persons in group A move to group I with a rate constant of γ.

Persons in both groups A and I move to the “Removed” group R at the same rate δ. This assumption is subject to debate. It is observed that almost all asymptomatic COVID-19 patients recover. In contrast, a fraction of infected patients die, and the rest recover. In this context, it is important to recall that the symbol R does notdenote “recovered,” but “removed,” either through recovery or through death. The point is that persons in group R are no longer susceptible. Therefore it is reasonable to assume that the rate at which persons in the I group move to group R and are thus “removed,” either through recovery or death, is the same as the rate at which persons in the A group recover and thus move to the group R. With reasoning, we use the same rate constant δ to multiply A and I.

These assumptions allow us to derive closed-form solutions to the simplified SAIR model, analogous to Theorem 1 for the SIR model. Define a new variable . Then it readily follows from (32) thatwhich is just the “SIR” model of (1) with M playing the role of I. Therefore Theorem 1 applies with M replacing I. Define, as before, the basic reproduction ratio σ ≔ β/γ. Then, in analogy with (5), with we have thatNow it follows from Theorem 1 that if σM 0 ≤ 1, then M(t)↓0, while if σM 0 > 1, then M(t) initially rises before decreasing to zero. In this case, the limiting value S ∞ can be computed from (3). Note that, because it follows that both A(t) and I(t) approach 0 as t → ∞. This is consistent with Theorem 6, which holds even without the simplifying assumptions (31). However, while (34) can assist in computing we would still wish to compute A and I individually. We now show this can be done.

One can substitute from (34) into (33) and observe that to get a differential equation that involves only S, namelywhere the term inside the square brackets is just some constant. Once S( · ) is found from (35), one can find M( · ) by simply substituting into (34). Next, upon noting that one can infer from (32) that

SAIR Model with vital dynamics

Until now we have studied the SAIR model without vital dynamics. Next we incorporate vital dynamics into the simplified SAIR model, as follows: Fig. 8 contains a flowchart of the above set of equations.

Fig. 8

Flowchart of the simplified SAIR model without vital dynamics.

In these equations, β is the infection ratio, α is the birth and death rate, γ is the rate at which group A and group I move to group R, and δ is the rate at which group A moves to group I. If we define and gather constants, we getThis is the same as (11) with I replaced by M.

The introduction of vital dynamics leads to a nonendemic equilibrium where the total infected population equals zero, and under suitable conditions, an endemic equilibrium. The nonendemic equilibrium is . As for the endemic equilibrium, S* and M* can be determined as in (12) and (14) with I* replaced by M*; thusThus an endemic equilibrium exists only when and the effective reproduction ratio is . Once M* is determined, we can compute A*, I* by setting which gives or

Now we establish the global asymptotic stability of the SAIR solution. The conclusion is analogous to Theorem 4. However, the method of proof is entirely different.

Define the set Then, whenever the trajectory (S(t), A(t), I(t)) of the system (37) converges to the endemic equilibrium (S*, A*, I*).

In analogy with (17), define the Lyapunov function candidateViewed as a function of (S, A, I), this function is convex because V is convex in (S, M) and M is linear in (A, I). However, V is positive semidefinite because it has its global minima along the line and not at the single point (S*, A*, I*). Next, in analogy with earlier arguments, it follows thatThereforeNow let us see what trajectories of (37) lie in the set . If for all t, then for all t ≥ 0. In turn this implies that orTherefore the only trajectories of (37) that lie in the set have for all t. It now follows from Theorem 5 that S(t) → S* and as t → ∞. Next, let us rewrite the equation for asHence, if M(t) → M* as t → ∞, it is a standard exercise in linear system theory to show thatThen it is a ready consequence that as t → ∞. □

In Robinson and Stilianakis (2013), the authors introduce vital dynamics into the SAIR model (26) without the simplifying assumptions (31). They derive formulas for the endemic equilibrium and establish that it is locally asymptotically stable by linearizing the model around this equilibrium. As of now, the problem of introducing vital dynamics into (26) and estabishing global asymptotic stability is still open.

We conclude this section by presenting a formula for the “basic reproduction ratio” for the SAIR model. In the SIR model with vital dynamics as in (12), the quantity σ defined in (13) is significant in that if σ < 1, then endemic equilibrium does not exist. Similarly, in the SEIR model with vital dynamics, it follows from (20) that an endemic equilibrium does not exist unless

***

In van den Driessche and Watmough (2008, Section 6.4), a very general solution is given for the basic reproduction ratio for a wide variety of models, which includes both the SEIR and SAIR models. Specifically, in van den Driessche and Watmough (2008, Section 6.4.2), the authors study the SEIR model with “small” secondary infections, which can be interpreted as an SAIR model. Therefore one can replace the term ϵβES in van den Driessche and Watmough (2008, Section 6.4.2) by and the results still hold. Specifically, after adjusting for current notation (and noting that in that paper), (van den Driessche & Watmough, 2008, Eq. (6.6)) becomesfor the SEIR model, while (van den Driessche & Watmough, 2008, Eq. (6.6)) becomesfor the SAIR model. In each case, the basic reproduction ratio σ must exceed one in order for the pandemic to increase initially, before subsiding.

Compartmental SAIR models

In large and diverse societies, it is not realistic to model the entire society as one homogeneous unit. It makes more sense to divide the society into a set of relatively homogeneous regions, which we refer to as compartments, and create models for each. In such a situation, the possibility of migration from one region to another is a distinct possibility, whatever be the “lockdown” policies in effect. One possibility is to divide the entire country into m compartments, and to create an overarching model. This would lead to a model with an enormous number of parameters to be estimated. Instead, we adopt what might be called the “thermodynamics” approach, wherein each compartment is deemed to interact with the rest of the country, often referred to as the “universe.” For another approach to the problem of migration, see Kaushal et al. (2020).

The two-compartment “thermodynamics” model without vital dynamics is as follows: This is just the general (not simplified) SAIR model of (26) with two extra terms: It is assumed that there is a migration from the universe U to the main compartment with a migration rate of μ, and similarly, there is a migration from the main compartment to the universe with a rate of μ. If both migration rates are zero, then we get two isolated SAIR models. Note too that migration is permitted only from the A and A groups. Clearly no country would permit migration from the infected groups. It is possible to make the above model more complex by permitting migration also from the S and S groups. It is left to the reader to show that Theorem 8 readily extends to this case as well.

Define the set in analogy with (28) , and define Then the trajectory of (42) and (43) approaches as t → ∞.

In analogy with Theorem 6, define the Lyapunov function candidatewhich is positive definite on . Thenwhich vanishes only on the set . Moreover, consists of the set of equilibria of the coupled system, and is thus an invariant set. The desired conclusion now follows from Theorem 7. □

Parameter estimation

In this section we revisit the SAIR model without vital dynamics (that is, without natural births and deaths), and show how the various parameters can be estimated.

Recall the general SAIR model:Note that the only variables we can observe are I and R. Since R stands for “removed” and not “recovered,” we could express R as the sum where H denotes the fraction that recover, and D denotes the fraction that die. Both quantities can be measured separately. From these observations, we aim to estimate the various quantities in the model.

In order to simplify the problem of estimating the parameters, we make a few assumptions.

It is assumed that . In other words, it is assumed that contact between persons of the S and A has the same likelihood of leading to infection as contact between persons of the S and I groups. There is some evidence to suggest that indeed viral shedding by asymptomatic persons is pretty much the same as that by infected persons; see for example Wölfel et al. (2020), He et al. (2020), Liu et al. (2020), Li et al. (2020). There also do not appear to be any strong biological arguments to say why this should not be the case.

In the case of asymptomatic patients, practically all of them recover, and very few if any die. In contrast, most infected patients recover, but some die. One could try to capture this situation by writingIn this notation, it will certainly be the case that γ > γ , or to put it in words, the fraction of A who recover is higher than the fraction of I that recover. However, it is assumed here thatthat, is the recovery rate for the A group is the same as the removal rate for the I group, which is the sum of the recovery rate and the death rate. With this assumption, we can write

Strictly speaking, in the model (44) we should incorporate a delay term to reflect the incubation period of the infection. However, as shown in Anderson and May (1991), Keeling and Rohani (2008) and reiterated in Robinson and Stilianakis (2013), the introduction of a delay term “rarely results in qualitatively different dynamics.” In Section 6, it will be seen that in the case of France and Switzerland, the lack of a delay term leads to worse estimates. However, this appears to be the exception rather than the rule.

With these simplifications, we arrive at the simplified SAIR model, namelyin addition to (45). It can be seen that there are only three parameters to be estimated here, namely β, γ, δ. Next we discuss the parameters can be estimated, in the order γ, δ, β.

To estimate γ, we start with the second part of (45), namely . Both R and I can be measured. However, some care is needed in estimating γ. If we were to estimate at several time instants (say by first-order differences), use I at the same time instants, and then use some kind of least-squares fit for γ, the results would not be very good. The reason is that in reality R is updated at discrete instants in time (usually once a day), and on top of that, the numbers can be adjusted up or down due to “reconciliation” of data. Instead, it would be better to write this relationship in integral form, asIn effect one makes a “phase portrait” of R versus I. Computing γ using a least-squares approach with the above relationship for various values of T gives an estimate that is more robust to the discrete nature of R.

Next we derive a method to estimate δ, using the data after lockdown. Suppose a “perfect” lockdown is implemented at time T, which causes after that time. In this case, the simplified SAIR model becomesIn particular, it follows thatAlso, the last equation in the model can be rewritten asSubstituting the expression for A( · ) and solving for I( · ) givesTherefore, ideally the plot of should be a straight-line with intercept log (δA(T)) and slope . By computing the slope we can estimate and by combining this with the earlier estimate for γ, we can get an estimate for δ. Note in passing that, once there is an estimate, it is possible to determine the fraction of asymptomatic patients by setting

Finally we come to estimating β, which turns out to be the most involved part. The solution methodology is based on the “closed-form solution” approach already introduced in Section 3. This system of equations describing the simplified SAIR model can be solved for pre-lockdown situation in terms of the reproduction rate by defining and observing that before lockdown, we havewhich can be solved in terms of aswhere and R(0) is assumed to be 0. Substituting the expression for M from (47) into (44) gives us the parametric solution in implicit form asand for the corresponding S, the infections M and I are found using :If the logarithm is approximated as then we get an early time solution for S in explicit form asandIn the early phase, we assume that A and I are approximately equal, so that M ≈ 2I. This gives the following expression for the early time infections.:

Since we can measure I(t) as a function of t, and we have an estimate for γ is available at this point, the above equation can be used to estimate β.

Numerical results

Now we present the outcomes of applying the parameter estimation techniques. In the first section, we apply our methods to eight countries from around the world, while in the second section, we analyze the situation in Delhi, which is (unfortunately) emerging as a “hotspot” in India.

Analysis of eight countries

Table 1 shows the estimated (with units of days) for various countries.

Table 1

Value of removal (recovery+death) time for different countries .

Country	γ^−1	Country	γ^−1
USA	50 ± 3	Brazil	20 ± 1
Italy	30 ± 2	India	20 ± 2
Iran	11 ± 1	Japan	11 ± 1
France	21 ± 2 or 100 ± 4	Switzerland
35 ± 2 or 19 ± 1

Fig. 9 shows the values of γ for various countries, by plotting R(T) versus as a function of T. It can be seen that, for six out eight countries, the plot is nearly linear, thus indicating a robust estimate for γ. However, for France and Switzerland, the graph is far from linear. We believe that this is because we ignored the incubation period of the virus. If we were to take this into account, then we would have to modify the dynamics as a delay-differential equation, in the formwhere τ is the incubation period. Thus

Fig. 9

Removal (recovery+death) frequency γ for various countries.

Fig. 10 shows the estimates of γ for these two countries using various values for the delay τ. It can be seen that, as τ is increased, the plot becomes more linear. It is not clear why this should be an issue only for two out of eight countries.

Fig. 10

Delay difference plots for : France(left) and Switzerland(right), in order to find the correct γ .

Fig. 11 shows the outcomes for estimating δ for various countries.

Fig. 11

Time is 14th February 2020. Extracting δ (using the decaying exponential in the analytical solution 10–30 days post lockdown).

Fig. 12 shows the outcomes for estimating β for various countries. This figure shows clearly that the “lockdown” has been implemented with quite varied levels of thoroughness in different countries.

Fig. 12

Time is 14th February 2020. Extracting β (using the analytical solution at early time).

Once we have fitted the parameters, we have solved the simplified SAIR model to generate the trajectory of the pandemic. Thus we have fitted the past history and made future forecasts for various countries in Fig. 13 . The quality of our estimates can be seen in this figure.

Fig. 13

Time is 14th February 2020. Numerical solution to the SAIR model using the estimated paramters β, δ and γ.

Table 2 presents the parameters of our models for all countries, in a convenient form.

Table 2

Parameters extracted by fitting the analytical solutions to the model we developed to the 7-day average data from the different countries. The ”–” (blanks) indicate that the parameter could not be estimated as country has not yet entered the region of exponential decay post lockdown.

Country	β	γ	δ	σ
USA	0.250 ± 0.02	0.020 ± 0.001	–	12.5 ± 1.625
India	0.242 ± 0.03	0.0486 ± 0.002	–	4.979 ± 0.821
Brazil	0.351 ± 0.03	0.0484 ± 0.003	–	7.252 ± 1.066
Iran	0.52 ± 0.02	0.086 ± 0.004	0.031 ± 0.001	6.046 ± 0.338
France	0.280 ± 0.012	0.046 ± 0.004	0.010 ± 0.001	6.08 ± 0.789
Switzerland	0.290 ± 0.009	0.0502 ± 0.001	0.071 ± 0.002	5.77 ± 0.321
Italy	0.256 ± 0.021	0.034 ± 0.002	0.091 ± 0.007	7.529 ± 1.054
Japan	0.185 ± 0.013	0.086 ± 0.003	0.022 ± 0.001	2.151 ± 0.225

Analysis of Delhi

The same methods were applied to analyze the situation in Delhi. Fig. 14 shows the various plots. The following are the key conclusions:

Fig. 14

Estimates and Predictions for Delhi: (a) Estimating γ. (b) Estimating β. (c) Reconstruction of past trajectory. (d) Prediction of future trajectory.

The impact of the lockdown in reducing both β (from 0.26 to 0.09) and σ (from 8.1 to 2.8) can be inferred from the data.

Future predictions of the progress of the disease show that “herd immunity,” in terms of the number of active infections peaking, will be achieved when the total of asymptomatic and infected is around 25% of the population. Given that the estimated value of σ after the lockdown is 2.8, classical SIR theory, based on ignoring asymptomatic patients, would predict that herd immunity is achieved at a level of or 64.29%. This shows that the computation of herd immunity must be modified for the SAIR model, and that the level for the SIR model is overly pessimistic. However, as of now, there is no explicit expression for the onset of herd immunity in the case of the SAIR model.

The above prediction completely ignores any kind of advances in the treatment of the disease. Obviously, the predictions will turn out to be overly pessimistic if any advances are made in prevention and cure of COVID-19.

Discussion and future research

In this paper, we have attempted to achieve two objectives. In the first objective, we undertook the task of completely analyzing the SAIR model which was introduced in Robinson and Stilianakis (2013) to incorporate asymptomatic patients. As a part of this, we established the global attractivity of the equilibria in the SAIR model, both with and without vital dynamics. Further, we extended the SAIR model to a compartmental model to accommodate migration. The major difficulty with the SAIR model is that it is not possible to observe the asymptomatic patients. Therefore, we provided a method for estimating the various parameters in a simplified SAIR model. The second objective was to validate our model by fitting the observed data. Our analysis shows that the model built upon our estimated parameters does an excellent job of explaining the evolution of the pandemic across several countries. We have also applied a similar analysis to the situation in Delhi. This analysis shows clearly the impact of implementing the “lockdown” in the Delhi area.

There is no shortage of interesting open problems to be tackled. We list some of them below:

Nonlinear Observers for the SAIR model: Perhaps the most interesting one is to design a nonlinear observer for the SAIR model. In the SIR model, which is the simplest, the medical system can measure both Infected (I) and Removed (R) populations; since in effect S can also be inferred. Hence the SIR model corresponds to a system in which all states can be measured. In the SEIR model, it follows from the above logic that the sum can be inferred from measurements of I and R. There is not much incentive to infer the values of S and E individually, because in the SEIR model, it is assumed that contact between members of the S and E groups does not lead to fresh infections. The situation is different in the SAIR model. As before, it is possible to infer the value of from the measurements of I and R. But the difference now is that there is a requirement to infer S and A individually, because contact between these two groups does lead to fresh infections. In control theory, there is a well-developed method of designing observers for inferring state variables that cannot be measured directly. While the theory is quite complete for linear systems, there are some results for nonlinear observers as well. It would be worthwhile to develop such observers for the SAIR model.

Lyapunov stability analysis of SAIR and compartmental models with vital dynamics: Several methods of stability analysis in the literature are based on “closed-form solutions” of the system at hand. The Lyapunov functions proposed here in Theorems 6 and 8 are new and can be applied to arbitrarily large concatenations of systems. However, this class of Lyapunov functions cannot be readily extended when vital dynamics are present. This problem is worth studying, because the introduction of vital dynamics actually makes the more realistic, by eliminating a continuum of equilibria, and resulting in just a few isolated equilibria.

Refinements of the SAIR model: It is possible to develop still finer models of the pandemic by introducing additional categories such as Quarantined, Healed, Ailing, Recognized (or Detected), Threatened, etc. The paper (Park, Cornforth, Dushoff, & Weitz, 2020) proposes what might be called an SEAIR model, but the level of analysis is not nearly so thorough as it is in the present paper. In Giordano, Blanchini, and Bruno (2020), eight different categories are introduced. By introducing more categories, we will get a more realistic model of disease progression. On the other hand, the number of parameters to be estimated increases drastically. The ideal trade-off between these two conflicting considerations remains to be explored.

Sensitivity to estimation errors: Given that the estimates for the various parameters are based on rather noisy and unreliable data, it would be desirable to carry out simulation studies to ascertain how sensitive the conclusions are to these error sources. For instance, in the estimates for Delhi, the values of σ (the basic reproduction ratio) are well away from 1 both before and after the lockdown. Hence we can be sure that the lockdown has had a beneficial effect, and moreover, this conclusion is robust against errors in the data and the consequent errors in parameter estimation. When it comes to applying epidemic models to real data, the literature is in the starting phase. We cite (Lavezzo, Franchin, & Ciavarella, et al., Verity, Okell, Dorigatti, et al., 2020) are two early examples. Another paper (Day, 2020) states that the COVID-19 was eliminated in a small Italian village thanks to identifying and isolating asymptomatic patients. The same paper also has the highest estimate currently available (75%) for the fraction of asymptomatic patients. For the most part, such applications are based on SIR models. As more and more such papers appear, it would be worthwhile to quantify the improvements in the accuracy of the predictions made by using SAIR-like models instead of SEIR-like models.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.