key: cord-0151139-z3drjvvl
authors: Borkar, Vivek S.; Manjunath, D.
title: Revisiting SIR in the age of COVID-19: Explicit Solutions and Control Problems
date: 2020-10-13
journal: nan
DOI: nan
sha: 278fa0c2025e05d64eadd1d947e4a8ec0803fc80
doc_id: 151139
cord_uid: z3drjvvl

The non-population conserving SIR (SIR-NC) model to describe the spread of infections in a community is proposed and studied. Unlike the standard SIR model, SIR-NC does not assume population conservation. Although similar in form to the standard SIR, SIR-NC admits a closed form solution while allowing us to model mortality, and also provides different, and arguably a more realistic, interpretation of the model parameters. Numerical comparisons of this SIR-NC model with the standard, population conserving, SIR model are provided. Extensions to include imported infections, interacting communities, and models that include births and deaths are presented and analyzed. Several numerical examples are also presented to illustrate these models. Two control problems for the SIR-NC epidemic model are presented. First we consider the continuous time model predictive control in which the cost function variables correspond to the levels of lockdown, the level of testing and quarantine, and the number of infections. We also include a switching cost for moving between lockdown levels. A discrete time version that is more amenable to computation is then presented along with numerical illustrations. We then consider a multi-objective and multi-community control where we can define multiple cost functions on the different communities and obtain the minimum cost control to keep the value function corresponding to these control objectives below a prescribed threshold.

Epidemic models have been used in a variety of situations-modeling of the spread of an infection in a community, spread of rumors in a community, spread of viruses in systems of computer networks, etc. An important workhorse in this domain has been the SIR model that was introduced by Kermack and McKendrick [17] and has received wide textbook treatment, e.g., [7] . This classical model assumes a community where everyone interacts with everyone in a homogeneous manner and the average behavior of these interactions are known. In contrast, in some new applications, the community interactions are represented by a graph and the interactions happen along the edges of the graph. See for example [10] for models of SIR epidemics on graphs. There are also stochastic models that explictly account for randomness, see, e.g., [2] , [3] . The standard SIR model for a single community is a compartmental model in which the population is divided into sub-populations of susceptible (S), infected (I), and recovered (R) individuals. At time t, let S(t) be the size of the population that is susceptible to be infected by the epidemic, I(t) be the number in the population that is currently infected and carrying the infection, and R(t) be the number that have recovered from the infection and are neither going to infect others nor will be infected again. Let µ be the contact rate 

Parameter λ is also called the pairwise infection rate. Observe that in this model, deaths, hospitalizations, and isolations are not modeled. In fact, it is assumed that those that are infected will continue to spread the infection during the entire time that they are carrying the infection, i.e., till they move into the recovered state.

Since the introduction of the SIR model, a myriad variations that include many more compartments and more detailed interactions have been proposed and studied. A key bottleneck in the study of these models is the nonlinearity that defies easy solutions. Till recently, explicit solutions of the basic system defined by (1) were not available. For example, [7] only provides a rather complex approximation for R(t). Solutions for some special cases are available, see, e.g., [22] . More recently, a parametric solution for the general case has been described in [12] . Even here, converting the parametric form to a time domain solution would require numerically solving an integral equation. Furthermore, small changes to the model that would make the model more realistic makes the search for a solution even more hopeless, as was exasparatingly noted in [19] 1 .

In the absence of explicit solutions, the focus has been on numerical evaluation and also on some key quantities of interest. For example, in the SIR model of (1), it is easy to see that λS(0) > γ will lead to a breakout of the epidemic. It is also known that at the end of the epidemic, the number that will not be infected will be strictly greater than 0, i.e., S(t) → c > 0. Several other key parameters are also of interest, especially if the infections result in the need for hospitalizations and other healthcare resources. The popular phrase in the context of the COVID-19 pandemic has been "flattening the curve", i.e., to delay and/or reduce the peak I(t). Using the parametric solution of [12] , we can get a closed form expression for the maximum value of I(t) (see (7) in a later section) in the SIR system of (1). However, the time at which this peak is achieved has to be obtained through a numerical solution to an integral equation.

Furthermore, in the absence of explicit solutions of the system of equations or of the key performance indicators, control problems are even harder. Thus, not surprisingly, application of interventions by suitably modulating either of λ or γ to achieve specific control objectives is not well studied. 2 In this paper we first make an important and, what we believe to be, a more realistic change to the basic SIR model in which we drop the population conservation assumption of (1) and propose the SIR-NC (for 'Non-Conservative') epidemic model. This is motivated by first recognizing that those that recover neither infect the others nor can they be infected by virtue of having developed immunity. (The SIR model confers immunity on those that have recovered. As of the writing of this paper, 1 While explaining that removing the population conservation constraint would make solutions for the even simpler SIS model impossible, the authors remark "It would seem that a fatal disease which this models is also not good for mathematics". 2 We remark here that this kind of control is different from vaccination based control whose effect is to reduce the susceptible population. This has received significant attention in the research literature. reinfections in COVID-19 are rare and the SIR is a very good model. The SIS model and its variants may be used for cases where re-infections are more common.) Thus we can treat R(t) as the part of the population that is 'dropping out' or removed from the system.

The outflow from the infected population could also include those that have died due to the infection, those that are cured and are free from the infection, and also those that are identified to be infected and hospitalized and/or quarantined. The parameter γ would then be the sum rate of all these. If one wants to focus on, say, the recovered, denoted by R(t), one would then haveṘ(t) = βI(t) where β < γ is the recovery rate. Similar formalism will be posited for the dead, the quarantined, etc. Thus our SIR-NC model allows us to model non-pharmacological interventions (NPIs) such as early isolation and quarantining of those that are infected through testing. Importantly, this is achieved without the addition of any new compartments to the model. The rest of the paper is organized as follows. In the next subsection we introduce the SIR-NC model and present its solution in the following section. In Section 2.1 we develop the model and present Theorem 1 that provides an explicit expression for (I(t), S(t)), the peak value of I(t), and the time at which the peak is realized. The case of time-varying λ and γ are also included. Theorem 1 is proved in Section 2.2. Several variations of the basic model are discussed in Section 3-analysis of a community with imported infections (Section 3.1), interacting communities (Section 3.2), and the SIR-NC model that allows for natural births and deaths (Section 3. 3) . Numerical examples will be used to illustrate these models. In Section 4 we discuss several control problems over a finite horizon for the SIR-NC model after discussing the control objectives, specifically, minimum cost control problems that control I(t) so as to achieve specific objectives. Sections 4.2 and 4.3 develop the finite horizon continuous time and discrete time control problems. In Section 4.4, we develop a more general multi-objective, multi-community control problem. Illustrative numerical results are also presented.

We begin by recalling the philosophy behind the classical SIR model. This model assumes that the population is comprised of three sub-populations of susceptible, infected, and recovered people. Let N (t) = S(t) + I(t) + R(t) denote the total population at time t.

There are no births or deaths, based on the assumption that these are on a much slower time scale and can therefore be ignored. When a susceptible person meets another person, nothing happens unless the latter is infected; the fraction of such meetings with an infected individual is I(t) N (t) . The infection happens with a rate of λ > 0, whence the rate of transfer from susceptible to infected is λ S(t)I(t) N (t) . The infected get recovered at rate γ > 0. The combined dynamics thus isṠ

Summing the three yieldsṄ (t) = 0, whence N (t) ≡ a constant N (say). The above then simplifies to (1) . While this classical model has been an extensively studied and successful model in epidemiology ( [13] has an excellent review, [7, 16] are examples of excellent textbook treatments) some of its limitations stand out: 3. Even in its pristine version, the dynamics is not easy to analyze and is usually studied only numerically or, as we have said previously, using messy approximations and hard to use parametric forms.

Motivated by this, we present here an alternative model that goes back to [11] and offers two clear advantages: implicit in the SIR model above, viz., that the recovered neither infect nor get infected (i.e., they become immune). They affect the evolution of S(t) and I(t) only through the probabilities of pairwise meetings between susceptible and infected, because R(t) enters the normalizing factor N (t) ≡ N . In our model we count only those pairwise meetings that matter, i.e., the meetings wherein a susceptible agent meets another susceptible or infected person. This is legitimate because the constraint that the total population is conserved is now dropped, whence the recovered can be legitimately viewed as having 'dropped out'.

We add the additional proviso that a susceptible person meets a large number of susceptible or infected persons-not necessarily the whole population, but a large enough statistically representative portion thereof. Then the probability of a susceptible person meeting an infected person is I(t)

where N (t) := S(t) + I(t). By a reasoning analogous to the above, the equations now become (with an additional tweak we mention later)

The additional change we have made is to replace γ by 0 < β < γ in the third equation.

That is, we allow for deaths due to the disease (more generally, including quarantine etc.)

at rate γ − β. Clearly, N (t) cannot be expected to be a constant, because the right hand side does not add to zero. This is not a complication, but a blessing in disguise, because now R(t) does not affect the evolution of S(t) and I(t) at all, in conformity with our view that the recovered 'fall out of the system'. Furthermore, the interpretation of λ in this model is different. This is described below.

Even before we discuss this any further, we wish to underscore the fact that this is not a case of sub-sampling random pairwise meetings that happen only one at a time on the microscopic scale in the standard SIR model. That would lead to a different dynamics. In other words, suppose that the individual meetings happen on a common Poisson clock with a fixed rate and we choose to observe only those between the susceptible and the infected or the susceptible. Any reasonable continuum limit of this should also re-scale γ as γ/ N (t) because of the time scaling implicit in the above sampling. Our model is distinct from this scenario in that it requires each person to meet a statistically significant portion of the entire population 'at one go'. This is, for example, what would happen at the workplace, in the school, during the commute, etc. This is also the view in the agent based simulation model implemented in [15] . What this stands for is that the 'atomic' or individual pairwise meetings take place on a faster time scale than the death rate due to the epidemic, which is a reasonable assumption. A simple scenario would be the discrete dynamics S n+1 = S n − ǫS n ζ n ,

where ǫ is the discrete time step, ζ n is an indicator random variable corresponding to a pairwise meeting which is 0 if it is benign (i.e., susceptible meeting susceptible) and 1 if not (i.e., susceptible meeting infected). The conditional probability that ζ n = 1 given the past is proportional to the fraction of infected population, say λ In Sn+In . Thinking of ǫ as a small time step, γ remains removal (or recovery) rate per unit time as ǫ ↓ 0, but the random meetings {ζ n } start crowding on the scale Θ 1 ǫ per unit time. Then the averaging effect kicks in and we get the desired dynamics in the limit. To make this precise, define

with linear interpolation on [nǫ, (n + 1)ǫ] for n ≥ 0. Then we have:

This is a special case of Lemma 1, p. 103, of [4] . In fact, this is the standard averaging effect of two time scale dynamics. Note that the interpretation of parameter λ gets modified from that for the original SIR model.

Getting back to our model, the most significant observation is what we have already made, viz., that S(t), I(t) obey a self-contained two dimensional differential equation which, as it turns out, is explicitly solvable as we shall see later.

2 Analysis of the SIR-NC model

In this section we analyze the SIR-NC model in detail. As we elaborated in the previous section, the R(t) only have a one-way coupling with I(t), the infected population. We thus have that

which can be explicitly computed once I(·) is known. Therefore the SIR-NC system that we need to analyze isṠ (t) = −λ I(t)S(t) S(t) + I(t) ,

Of course, we should have I(0) > 0 and define S(0) + I(0) = N (0). We could also have R(0) > 0, but since they do not affect the evolution, R(0) = 0 is without loss of generality.

This change to non-conservation allows us to provide explicit expressions for S(t), I(t), I max := the peak value of I(t), and T max := the time at which it is attained (Theorem 1).

This in turn helps us define several control problems with the many different control objectives to shape the evolution of the epidemic.

It is easy to see that the condition for an epidemic to break out in SIR-NC is

The main result in this section is the following theorem. (See also [11] .).

, the solution to the SIR-NC system is as follows.

For the more general case of time-dependent λ(t), γ(t), t ≥ 0, we similarly have 

Remark The last equation is arrived at simply by putting the time derivative of I(·) equal to zero. This need not have a unique solution, nor need every solution correspond to a maximum, even a local maximum.

We now compare the computations from this model with that of the standard popula-tion conserving SIR model. First, R(t) for the SIR-NC is obtained by using I(t) from (4) in (2) . Compare this with with the following approximation from [7] . Here ρ = λ/γ.

Next, compare the computation of T max in SIR-NC as above with that in the SIR model using the the parametric form of S(t) and I(t) from [12] . The following is adapted from [24] .

As can be seen, t does not have a closed form solution.

Finally, we consider the computation of I max . In the SIR model we can see that the maximum of I(t) occurs for u = S(0) λ γ in (6) and the I max in the SIR system can be evaluated as

where ρ = γ/λ. In the approximation above, we have assumed I(0) is very small compared to S(0) and hence S(0) ≈ N. The corresponding expression for the SIR-NC system is

where, as before, ρ = γ/λ and the approximation assumes that I(0) is very small compared to S(0).

We now discuss how the key differences between the SIR and the SIR-NC models affect 000  135  138  54  53  68  69  10,000  181  184  70  69  91  92  100,000  228  230  85  84  114  115   Table 1 : Table shows T max different combinations of N, λ, and γ in the SIR and SIR-NC systems. In all of these we assume I(0) = 1 and R(0) = 0. for the two models. It is also instructive to compare I max and T max for the two models;

Fig. 2 compares the Imax N for the SIR and the SIR-NC models as a function of ρ = γ/λ. As expected, I max in SIR-NC is higher than that in SIR. Table 1 lists T max for the SIR and SIR-NC systems for different values of N, λ, and γ. Interestingly, we see that T max in both the systems are very nearly the same. We performed these computations for several more parameters and this observation appears to be more generally true.

In the rest of the paper, we will only consider the SIR-NC model.

In Fig. 3 we show how I max and T max behave if the infected could be removed at a faster rate than the natural removal rate. The marginal relative increase in T max is significantly lower than the marginal decrease in I max for small γ 1 . In the example, observe that doubling be obtained by using (4).

We now consider more general cases that may be applicable in epidemic conditions.

We argue later that the following explicit form of λ(t) would be of interest.

where 0 < a ≤ 1. This models a limited lockdown period during which there would be a significant reduction in the contact rate and hence the infectivity. Further, we will use γ(t) = γ 0 (1 + bt) which models a super-exponential increase in testing capacity and isolation. In Fig. 4a , γ(t) is a constant (γ 0 ) and λ(t) has the form above. We plot I(t) for different values of t 1 . We see that an early lockdown (corresponding to t 1 = 15 and 25) only delays the T max but does not change I max . If the lockdown is started close to the peak, it results in 'second wave' with the second peak being lower than the first. A moderate delay in the start of the lockdown (t 1 = 40) arrests the first growth but does not significantly affect I max or T max . If the lockdown starts close to the peak, I(t) comes down sharply but causes a second peak! In Fig. 4b we plot I(t) keeping λ(t) constant, i.e., no lockdowns, and use γ(t) as above. This is to be expected because, the lockdowns only postpone the inevitable. bserve that even a small amount of increased testing and quarantining can be very effective in reducing I max significantly; a 1% increase in the rate of removal per unit time (doubling testing capacity in 70 time units) reduces I max by about 30% whereas a 3% (double capacity in 24 days) can decrease I max by 70%. T max is not affected significantly in this. Finally, in Fig. 4c , we plot I(t) when there is lockdown and also the γ is increased.

We see that a lockdown and increased testing can bring the epidemic under control.

The preceding motivates us to propose formal optimal control frameworks. They are detailed in Section 4.

Letting

,

we have x(t) + y(t) ≡ 1. Also, x(0) < 1 (by the assumption that I(0) > 0). We can rewritė

x(t) andẏ(t) as follows.

Thenẋ(t) +ẏ(t) = d dt (x(t) + y(t)) = 0 as expected. Noẇ

For γ < λ, x(t) → 0 and y(t) → 1 as expected, i.e., for x(0) ∈ [0, 1), (0, 1) is the asymptotically stable equilibrium. It also follows by setting t = 0 that

. We are now ready to derive the expression for I(t). Letting N (t) = S(t) + I(t), we havė

Then clearly

Then I(t) = y(t)N (t) → 0. Also, from the explicit expression for S(t), we have S(t) → 0, thus N (∞) = 0.

To determine T max , note that initially x(t) ≈ 1. Coupled with λ > γ, we haveİ(t) > 0, so that I(t) increases till it hits the isocline S N = 1 − γ λ at time T max . Thus T max can be calculated as follows. We havė

Thus T max is the solution of the equation

This completes the proof of the constant coefficient case. The proof for time-dependent λ, γ goes exactly the same way.

Now considerṘ

for some β < γ. Thus β is the recovery rate and γ − β the death rate. Then

We now consider several variations of the basic SIR-NC model that we introduced in the previous section. Specifically, we present two types of variations and characterize the solutions. The first set of variations is motivated by the actions of several communities during the COVID-19 pandemic. Communities are isolating themselves to different degrees and also applying differing degrees of social distancing rules within the community. To analyze this, in Section 3.1 we consider the case when a community has 'imported infections', i.e., an external agency introduces infections into the community at a steady rate. In Section 3.2

we use the results of Section 3.1 to analyze the interaction two separate communnities that interact with each other. We obtain explicit characterizations for both these models. The second variation is to introduce natural births and deaths in Section 3.3 where we analyze the equilibrium behaviour of the model.

In many communities, infections are also induced by interactions with other communities.

Assume that each of the susceptible in the population meets µ e persons from outside of the community and p e of those meetings result in an infection. Defining ν = p e µ e , we have the following SIR-NC model in which infections are also imported from other communities.

We observe that the condition for an epidemic to break out is

. Along the lines of Section 2.2 we can show that x(t) satisfies the

This can be written as

That is,

. Now,

We have just derived the following replacement for Theorem 1 for the case of SIR-NC with imported infections. T max for this model can be obtained as in Theorem 1.

Theorem 2. The solution to the SIR-NC system with imported infections is as follows.

where 

We now illustrate this model with some numerical examples. Fig. 5a illustrates I(t) for different values of ν. Observe that even a small ν advances T max significantly, ν = 0.01λ reduces T max by nearly half (from about 55 to about 29). This dependence is studied in more detail in Fig. 5b that shows T max and I max as a function of ν. The marginal reduction in T max is significant for small values of ν but I max is less dramatic.

Now consider two interacting communities, labeled a and b, that have different intracommunity contact and infection infection rates, and in addition, inter-community contact and infection rates. Let N a and N b be the size of the two communities. Defining λ a , λ b , γ a , γ b and λ ab in the obvious way, we have the following equations defining the system evolution.

Letting

we have the coupled dynamicṡ

and an analogous dynamics for (y a (·), y b (·)). This is a cooperative dynamics in the sense of [14] and by virtue of having bounded trajectories, will converge for almost all initial conditions to the set of its equilibria. These can be found by setting the right hand side equal to zero, and their stability can be analyzed using local linearization. We do not pursue this here. This equation, however, is not analytically tractable. So we aim for an approximate solution by assuming that λ ab = ǫ, 0 < ǫ << 1 T , where T > 0 is the time horizon of interest. Let S 0 i (t), I 0 i (t), i = a, b, denote the decoupled dynamics corresponding to ǫ = 0. Then for general ǫ > 0, we can use the results of the preceding section to obtain an approximation for S a (t), I a (t), T a max (:= the corresponding peaking time) as follows. However, the two are coupled. We can decouple them by resorting to the approxima-

An exact expression for this error term can be given in terms of Alekseev's nonlinear variation of constants formula [1] , [5] as follows.

View the dynamics (11))- (13) , (17) of community a (a symmetric argument works for community b) as a perturbation of the uncoupled dynamicṡ

The linearization of (19)-(20) around the nominal trajectory (S(·), I(·), R(·)) is given as

Write this asq 

Then by Alekseev's nonlinear variation of constants formula [1] , [5] ,

This is an exact expression for the net perturbation of the solution, albeit implicit. The second term on the right is O(ǫ). Thus if we ignore the O(ǫ 2 ) terms, we can write

This gives an explicit expression for the net perturbation of the solution, albeit only approximate to within O(ǫ 2 ). Note, however, that this approximation will be valid only for a finite time interval, because the hypothesis that ε(t) = O(ǫ) may not continue to hold unless T is small. For computational purposes, the above procedure can be repeated over such small time intervals.

Remark: The case of λ ab = λ ba is also of interest, e.g., when one community practices the safety recommendations more strictly than the other community. This case can be analyzed analogously. Fig. 6 illustrates the effects of the intercommunity interaction rates, λ ab and λ ba , on I(t)

when different kinds of communities interact. In Fig. 6a we consider two similar communities with the same initial population and same λ. The trajectories for I(t) for both communities are, as is to be expected, identical. However, even a small amount of interaction can significantly affect the trajectory, leading to increased I max and a reduced T max .

This is in line with the results of the preceding subsection. In Fig. 6b , we see that when λ a = λ b and λ ab = λ ba , and the sizes are significantly different, the trajectories of I(t) are very similar in both the communities and they peak at nearly the same time. However, if N a = N b and λ a = λ b while λ ab = λ ba , then even a small amount of interactions matter as can be seen in the significant shift in both T max and I max for community a which has a lower λ. This is seen in Fig. 6c . Finally, in Fig. 6d we show I a (t) and I b (t) when there is asymmetry in the two communities both in terms of the size and the λ. It is more interesting to see I a (t) + I b (t) for this last case, as shown in Fig. 7 . Observe that if the interaction is small, the total I(t) has a comparatively 'fatter' middle, possibly with a double hump.

This disappears when the interactions increase.

The models for long duration epidemics, like the way Covid-19 is expected to be, should include natural births and deaths. The standard SIR model has been extended to include natural births and deaths. However, almost all the models continue to make the population conservation assumption-the natural death rate is equal to the birth rate and the total population is conserved. The SIR-NC model naturally allows for unconstrained birth and death rates. Thus in this section we describe the SIR-NC model that includes births and natural deaths with no constraints on the rates. However, if we introduce some restrictions on the parameters, not unlike those that are derived from the standard SIR model, we can provide additional insights to the properties of the evolution of the epidemic. Here, we make the natural assumption that all births are susceptible.

Consider the dynamics with a small birth rate and a small death rate for 'normal' deaths (as opposed to deaths due to the epidemic), given by:

where N ǫ (t) = S ǫ (t)+ I ǫ (t)+ R ǫ (t) and κ, υ i , ν i = O(ǫ), i = 1, 2, for some 0 < ǫ << 1. Note that we have now included the recovered, i.e., R ǫ (t) in the definition of the normalizing factor N ǫ (t) in view of the fact that it directly affects the evolution of S ǫ (·) given by (26) and can therefore no longer be viewed as having 'dropped out'. The υ i 's are birth rates and the ν i 's the rates of normal deaths for the infected and recovered populations respectively, and are positive. The parameter κ is the net rate of birth minus normal deaths for the susceptible population, and can be either positive or negative.

In Fig. 8 we show the trajectories of S(t), I(t), and R(t) for this model for three different cases: birth rate higher than natual death rate, balanced birth and natural death rates, and natural death rates higher than birth rates. We show both the short term and the long term trajectories. Observe that in all the cases, there are significant oscillations in S(t), I(t), and R(t). The long term trajectory of I(t) for the three cases are shown in more detail in Fig. 9 . We observe that I(t) decays very slowly and I(3000) has values 219, 102, and 50 for the three cases of, respectively, excess births, balanced, excess deaths.

The preceding observations can be made more precise with the following additional assumption:

This in particular includes the case of

corresponding to a balance of birth rate and the rate of natural deaths community-wise.

The former also includes the case of excess of natural death rates over the birth rates. 

there is nothing to prove, so suppose N > 0. Then every limiting trajectory (S ′ (·), I ′ (·), R ′ (·)) of (26)-(28) in the ω-limit set Ω of (26)-(28) must satisfy I ′ (·) ≡ 0. But then (28) forces R ′ (·) ≡ 0, so that S ′ (·) ≡ N . But (0, 0, N ) is an unstable equilibrium for N > 0, since a small perturbation that renders I ′ (0) > 0 will make the trajectory move away from it. The claim follows.

Therefore a necessary condition for non-extinction is that at least oneof the subcommunities should have a non-zero value for rate of birth minus the rate of normal deaths,

i.e., κ > 0, or v 1 > ν 1 , or v 2 > ν 2 . The most interesting aspect of this analysis is the case ( † †), which subsumes the conditions under which the classical SIR model is derived. Recall that in the classical SIR model, non-extinction is possible. For SIR-NC, however, even a small death rate due to the epidemic leads to extinction.

More generally, we can say the following. Proof. If (s, i, r) is an equilibrium of the above system, we perforce have r = β ν 2 i by setting the right hand side of (28) to zero and therefore the equilibrium equations obtained by setting the right hand side of (26)-(27) equal to zero reduce to

Then for any non-zero equilibrium, i > 0 and s = ζi/κ where ζ := γ + ν 1 − υ 1 − υ 2 β ν 2 . This fixes the value of λ in terms of the other parameters. Thus for all but one value of λ, the dynamics does not equilibrate.

We now discuss another key objective of this paper-to define associated control problems.

The following objectives are of obvious interest to the controller.

• Increasing T max : this allows the healthcare system to prepare better to treat the infected population, especially those that may develop severe symptoms.

• Decreasing I max : This reduces the preparation cost.

• Decreasing T e , the time to end the epidemic: This would mean that the system can get back to normalcy quicker.

• Achieve herd immunity by T 0 .

In addition one could also define more detailed objectives, e.g., define objectives similar to the preceding for specific communities that may or may not interact.

Two parameters are available for control; λ := the spreading rate and γ := the removal rate. Recall that in the original SIR model, λ is the product of µ := the rate at which people meet (or more accurately, are in physical proximity), and p := the probability of the proximity resulting in an infection. In SIR-NC, we have separated the relative time scales of the two processes, but the interpretation remains. While p is more likely to be a feature of the infection and the community, µ can be controlled through measures such as imposing social distancing and varying degrees of lockdown. Hence we let λ take finitely many discrete levels. Further, as we have argued earlier, the recovered population is viewed as removal. In practice, increasing the removal rate is effected by identifying the infected early, i.e., by increased testing and quarantining of those that test positive. We assume that the testing capacity can be continuously varied and thus allow γ ∈ a suitable interval of ℜ + .

Specifically, we assume that λ ∈ L := {ℓ 1 , · · · , ℓ M }, i.e., takes discrete values, with switching cost (ℓ ′ , ℓ ′′ ) → c(ℓ ′ , ℓ ′′ ) > 0 when ℓ ′ = ℓ ′′ and = 0 when ℓ ′ = ℓ ′′ . It is customary and sensible to assume that c satisfies the 'triangle inequality'

for distinct ℓ i , ℓ j , ℓ k ∈ L, in order to avoid some pathologies. On the other hand, γ can be continuously manipulated as a function of the current state. We first consider as the state process y(t) = I(t) N (t) ∈ [0, 1] introduced earlier. Later on we shall also consider the full state (S(t), I(t)). We let γ : [0, 1] → [0, Γ] for some Γ ≥ 0.

We consider a finite horizon control problem. There is a very good motivation for doing so. For scenarios such as this, Model Predictive Control (MPC) [18] is often the paradigm of choice, where at time t, one solves a finite horizon optimal control problem on time horizon [t, t + T ] and uses the optimal control choice for time t at time t. However, as time moves forward, to say t ′ > t, one again solves a new finite horizon problem on [t ′ , t ′ +T ] and the process is repeated for each time instant. The basis of this procedure is the solution of a finite horizon problem, so we have preferred it over other alternatives.

We begin with the continuous time formulation. Let T > 0 be the time horizon under consideration. There is a per unit time running cost k(y(t), λ(t), γ(t)) ≥ 0 and a terminal cost h(y(T )) ≥ 0. The objective is to minimize the cost T s k(y(t), λ(t), γ(t))dt +

for s = 0. Here {τ i } are the times when λ(t) makes a discrete switch between two values in L, λ(τ − i ), λ(τ i ) denoting resp., the value of λ(·) just before the switch and just after the switch. The value function then will be a map

defined as the infimum of (29) over all admissible choices of λ(t), γ(t), τ i , T ≥ τ i ; τ i , t ≥ s, with y(s) = y ′ , λ(s) = ℓ i . For ease of notation, we shall write V (y, ℓ i , t) simply as V i (y, t) and the corresponding running cost k(y, ℓ i , γ) as k i (y, γ). Since this is a mixed switching and continuous control problem, the Hamilton-Jacobi system here is a system of quasivariational inequalities given by: for all i,

The explanation is as follows. In absence of any switching of the λ parameter, (30) with equality would be the standard Hamilton-Jacobi equation of dynamic programming for continuous control. The same will be true now when you don't switch. When it is better to switch, continuing with the continuous control without switching will not be optimal, so inequality should hold. This explains (30). Inequality (31) is in similar spirit for the discrete control decision of switching: equality holds when it is optimal to switch, inequality holds otherwise. At any instant, it is either optimal to switch or to continue with the optimal continuous control. This leads to (32). Finally, (33) is the boundary condition given by the terminal cost.

Given a solution of this system, one uses λ(t) = ℓ i when

along with the continuous control γ(t) := a measurable selection from the minimizers of

Since we are in a finite horizon, a discretization can be solved by simple backward recursion.

Next we state an alternative framework based on the actual state variables, i.e., the pair (S(t), I(t)). The advantage of using y(t) as state variable is that it is a bounded scalar and therefore easier to handle. It also reflects the costs associated with I(t) somewhat faithfully. Nevertheless, it is not the real thing and computational power permitting, it may be better to stay with the actual variables. We shall use the same notation for the cost functions and the value function as before. The quasi-variational inequality then becomes

Comments analogous to those following (30)-(33) also apply here.

Unfortunately, the quasi-variational inequalities for first order systems are difficult to analyze and are usually not well-posed in the classical sense. For this reason, most work in this direction has been in the framework of viscosity solutions [9, 8, 25] . We shall take recourse to a simpler alternative, viz., to consider a discrete skeleton of the problem, which leads to a simple backward recursion.

We should mention here that there is already some work on pure impulse control of this model, i.e., without the continuous control [6] , [20] . This is technically a little easier and does not run into some of the aforementioned technical issues.

Let a > 0 be a small time-discretization step. A discretization of the underlying dynamics is given by the Euler scheme for the original continuous dynamics as:

for t = 0, 1, 2, · · · . We again use the same notation as before for cost functions and the value function. The objective is to minimize 

for s = 0. The value function V i (S ′ , I ′ , s) is then the infimum of (40) over all admissible λ(t), γ(t), τ i , T ≥ t, τ i ≥ s, beginning at λ(s) = ℓ i and S(s) = S ′ , I(s) = I ′ . The quasivariational inequalities can be written along the lines of the foregoing, but for discrete systems, they can be put in the form of a very convenient backward recursion. Let Z(t) = (S(t), I(t)) and write the above discrete dynamics (38)-(39) compactly as

Then the discrete quasi-variational inequalities reduce to the simple backward recursion of dynamic programming given by

with terminal condition V i (z, T ) = h(z).

Remark (40) can also be used for the population conserving SIR model with an appropriate version of the discrete dynamics for (41). We do not pursue this here.

If the objective is to minimize the peak over this interval, we augment the dynamics 

The new state now is Z(t) = (S(t), I(t), M (t)) with dynamics

for a suitably redefined F . The dynamic programming equation is the same as before with k i ≡ 0 ∀i and h((s, i, m)) = m.

We now illustrate the preceding control problem via numerical examples. In our computations, we minimize (40) subject to (41) with a time-discretization step of a = 1, treating it as an optimization problem, and use a simple forward tree search to solve it.

We consider a community with N = 10, 000, natural λ of 1/4, and natural γ, designated as γ 0 , of 1/15. Recall from the example in Section 2.1 that for this system T max = 55.4 and I max = 4523. We assume that the control objective at time t would be to have the number of infected at time (t + T ) to be a fraction α(t + T ) of the value in the uncontrolled system, which we will denote by I o (t + T ). With this view, we use the following cost function h(·) h(S(t + T ), I(t + T )) = A 1 e a 1 (I(t+T )−αIo(t+T )) , with a 1 A 2 to be constants. This function penalizes missing the the target significantly more than the reward for exceeding it. In our computations we have used a 1 = 10, A 1 = 100, and T = 3. To model an aggressive control objective we use α(t) = 0.1 for all t. To model a potentially more modest objective, we use

This corresponds to slowly increasing the requirement reaching 0.1I max at T max .

We assume M = 3, i.e., there are three contact rates among the members of the community corresponding to three levels of lockdown. We use L = {1/4, 3/16, 1/8} corresponding to, respectively, no lockdown, a partial lockdown, and a total lockdown. The switching costs are

The cost of switching to a lower λ is significantly higher than to switch to a higher λ. We use a 2 = 10, 000 in the simulations.

The running cost, k(·), will have the following additive components.

• Cost due to the chosen level of lockdown; we will use L 1 = 0, L 3 = 10L 2 corresponding to the costs of the three levels of the lockdown. We use L 2 = 10, 000.

• Cost of testing and quarantining. This will be an increasing function of the rate of removal of the infected subpopulation in excess of the natural rate of 1/15 and the number that are infected, I(t). At the beginning of the epidemic, it may be hard to obtain the testing kits and to also define the protocol and will hence be very costly to achieve a high rate of removal. Eventually, the time dependent cost becomes negligible. This leads us to propose the following function for this component: There is an interesting change in the behaviour of the cost functions and the control variable as a 3 is increased from 1.0. We will see this later.

From the preceding discussion, we have

In Figs. 10 and 11 we plot I(t), γ(t) and the running cost for different a 3 and α(t).

There are several takeaways from these results. For the parameters chosen, the switching costs and the running cost of reduced λ was too high for the system to switch λ, thus λ(t) = 1/4 was the control for all the parameters. As we changed a 3 from 1.0 to 1.2, for decreasing α(t), the change in both I(t) and γ(t) was barely noticeable. This was the case even when α(t) was held constant at 0.1. However as we changed a 3 from 1.2 to 1.35 the changes were significant both in shape and magnitude of I(t), γ(t), and the total running cost. This can be seen when comparing the plots in Figs. 10 and 11. A higher a 3 induces higher levels of control and hence reduced I(t) this in turn reduced the running cost. Observe that there is an order of magnitude reduction in I(t) and in the cost. Thus a higher 'valuation' for infections induces a higher level control which in turn reduces costs in the long run.

As we have seen at the beginning of this section, there are several competing criteria we may want to optimize and focusing on just one at the expense of others may not be wise. We therefore propose here a multi-objective control scheme inspired by [23] . We shall assume that the state space ℜ + × ℜ + is now suitably discretized and denoted by S. Likewise, the dynamics (41) is suitably replaced by a dynamics on S, where we retain the same notation as (41) for sake of simplicity. For example, the original F could be replaced by F followed by a suitable vector quantization scheme. As already noted, some cost criteria may require a suitable augmentation of the state space.

Consider K > 1 running cost functions k i 1 (·, ·, ·), · · · , k i K (·, ·, ·)), switching cost functions c 1 (·, ·), · · · , c K (·, ·) and terminal cost functions h 1 (·), · · · , h K (·). Let C 1 , · · · , C K > 0 denote prescribed numbers. Let (i, z, t) → V i j (z, t), 1 ≤ j ≤ K, denote the value function corresponding to the jth running cost function k · j (·, ·, ·), switching cost function c j (·, ·), and terminal cost h j (·). Our objective will be to ensure that for given initial condition z = z 0 of the state and the initial choice ℓ i 0 of the discrete control variable,

We assume that the set of admissible controls for which (44) holds is not empty. We scalarize this problem by considering a single running cost function (k · w (·, ·, ·) := j w j k i j (z, u, t), a single switching costĉa w (·, ·) := j w j c j (·, ·) and a single terminal cost functionĥ w (·) := j w j h j (·), where w 0 := [w(1), · · · , w(K)] T be a vector satisfying w(j) > 0 ∀j and j w(j) = 1 (i.e., a probability vector). The algorithm then is as follows. For n ≥ 0, do the following.

with terminal condition V i m,n (z, T ) = h m (z), 1 ≤ m ≤ K.

For n ≥ 0, do the following.

with terminal condition V i n (z, T ) =ĥ wn (z).

3. Perform a single update of w n = [w n (1), · · · , w n (K)] T as 

Denote by V i w the unique solution to

Iteration (48) is simply the Euler scheme for the differential equatioṅ

where V i (z, t) := m w m (t)V i m (z, t). Let (q n (·, ·), v n (·, ·))) : (ℜ + ) 2 × {0, 1, · · · , T } → L × [0, Γ], n ≥ 0, be a sequence of stationary policies such that (q n (z, m), v n (z, m)) attains the minimum on the right hand side of (47).

Theorem 5. Every stationary policy obtained as a limit point of (q n , v n ) as n ↑ ∞ satisfies (44). This is proved exactly as in [23] . In fact, the scenario here is vastly simpler than that of ibid. due to absence of noise, asynchrony, etc. We only sketch the key steps here, referring the reader to [23] for further details. First, observe that (50) is a special case the replicator dynamics of evolutionary biology [21] . As in [23] , the discrete iteration (which can be viewed as the Euler scheme for (50) with slowly decreasing stepsizes) tracks the asymptotic behavior of (50). The condition (49) plays a key role here: since a(n) plays the role of a time step in the time discretization of (50), (49) ensures that the entire time axis is covered, so that the asymptotic behavior of (50) can be tracked. We do not, however, need the standard condition in stochastic approximation algorithms that n a(n) 2 < ∞ which features in [23] , because we have a purely deterministic iteration without any noise.

Thus the problem is reduced to establishing that (50) indeed converges to the desired set. This is already proved in [23] .

Next we consider the control of two interacting communities like in Section 3.2 with the stipulation that the dynamics of one evolves on a much faster scale than the other. One possible scenario is a dense slum interacting with a geographically sparse affluent section.

Following Section 3.2, we model this as the coupled system

N a = S a (t) + I a (t) (55)

Here 0 < ε << 1 is a small parameter that renders the evolution of (S b (t), I b (t)) to be on a slower time scale. Following the standard philosophy for analyzing two time scale systems, we view the slow time scale dynamics as quasi-static for purposes of analyzing the dynamics on the fast time scale and treat the fast time scale dynamics as quasi-equilibrated for the purposes of analyzing the slow dynamics. Thus for the former, we consider the system dS a (t) dt = −S a (t) λ a (t) I a (t) N a (t) + λ ab (t)Z ,

where we view I b (t) N b (t) ≈ a constant Z. This coincides with our model of SIR-NC with imported infections. If we follow the standard procedure described above, then under the condition λ ab < λ a − γ a , the asymptotic equilibrium is the origin, i.e., extinction, which means that the fast variables simply drop out and the slow dynamics can be analyzed as a single community dynamics. To get something more informative than this, we propose the following, particularly in the context of the control problem. Choose ε < a << T above so that the time horizon T > 0 is divisible by a and a is divisible by ε. (If not, we can replace T /a, a/ε by ⌈T /a⌉, ⌈a/ε⌉ resp. in what follows.) Consider the discrete time intervals K n := {t(n), t(n) + εa, · · · , t(n + 1)} for 0 ≤ n < N T := T aε . On the interval K n , beginning with n = N T − 1 and going backward in time, do the following:

1. Consider the discretized dynamics for state variables (S a (n), I a (n)) corresponding to (57)-(58) for each possible value of Z := I b (t(n + 1))/N b (t(n + 1)).

2. Solve the dynamic programming equation on the discrete time interval J n := {t(n), t(n)+ a/ε, t(n) + 2a/ε, · · · , t(n + 1)} for the discretized dynamics of (57)-(58) with terminal cost := the time t(n + 1) value function inherited from the preceding computation on [t(n + 1), t(n + 2)] if t(n) < N T , and := h otherwise, so as to obtain the value function for m ∈ J n and initial (i.e., at time T (n)) condition (s, i) ∈ (ℜ + ) 2 .

Compute the optimal trajectory of (S a (m), I a (m)), m ∈ J n . 

We have proposed a variant of the classical SIR model wherein we drop the key assumption of conservation of population size. This Non-Conservative SIR or SIR-NC for short, has the advantage of being analytically tractable, leading to explicit expressions in the simple cases, and in addition, being more realistic for lethal pandemics. We also introduced several variations of the basic model, such as inclusion of natural birth and death rates, external inputs, interacting communities, etc. While some preliminary results have been presented for these, they open up many new directions for future work.

Finally, we consider controlled versions of the above. After introducing a variety of control objectives and arguing that Model Predictive Control is a natural framework for pandemic control, we propose a finite horizon control problem. First we do so for the original continuous time case, leading to a system of quasi-variational inequalities due to co-occurrence of continuous controls with discrete controls that also have switching costs.

We point out the technical difficulties therein and resort to a discrete version, leading to drastic simplification. We also consider multi-objective control in this framework. Finally, we briefly touch upon interacting communities with significantly separated time scales and propose an approximation scheme for analyzing the associated control problems. Once again, this opens up rich possibilities for future work.

An estimate for the perturbations of the solutions of ordinary differential equations

The threshold behaviour of epidemic models

A modification of the general stochastic epidemic motivated by AIDS modeling

Stochastic Approximation: A Dynamical Systems Viewpoint

Perturbations of non-linear systems of differential equations

Optimal intervention for epidemic models with general infection and removal rate functions

Epidemic Modelling: An Introduction

Hybrid control systems and viscosity solutions

Optimal switching for ordinary differential equations

The effect of network topology on the spread of epidemics

The spread of epidemics

Exact analytical solutions of the susceptible-infected-recovered (SIR) epidemic model and of the SIR model with equal death and birth rates

The mathematics of infectious diseases

Systems of differential equations that are competitive or cooperative ii: Convergence almost everywhere

Covid-19 City-Scale Simulation Team, Covid-10 epidemic: (unlocking the lockdown in india

Modeling Infectious Diseases in Humans and Animals

A contribution to the mathematical theory of epidemics

Model predictive control: Recent developments and future promise

An integrable SIS model

Population Games and Evolutionary Dynamics

A note on exact solution of SIR and SIS epidemic models

Q-learning for Markov decision processes with a satisfiability criterion

Susceptible-infected-recovered (SIR) dynamics of COVID-19 and economic impact

Optimal control of hybrid systems and a system of quasi-variational inequalities