Biology
Computer Science-2021.gz:143161
Computer Science
bc09f20c7820604f99b22123578c0d581aa2ca2c
SUTRA: A Novel Approach to Modelling Pandemics with Applications to COVID-19

The Covid-19 pandemic has two key properties: (i) asymptomatic cases (both detected and undetected) that can result in new infections, and (ii) time-varying characteristics due to new variants, Non-Pharmaceutical Interventions etc. We develop a model called SUTRA (Susceptible, Undetected though infected, Tested positive, and Removed Analysis) that takes into account both of these two key properties. While applying the model to a region, two parameters of the model can be learnt from the number of daily new cases found in the region. Using the learnt values of the parameters the model can predict the number of daily new cases so long as the learnt parameters do not change substantially. Whenever any of the two parameters changes due to the key property (ii) above, the SUTRA model can detect that the values of one or both of the parameters have changed. Further, the model has the capability to relearn the changed parameter values, and then use these to carry out the prediction of the trajectory of the pandemic for the region of concern. The SUTRA approach can be applied at various levels of granularity, from an entire country to a district, more specifically, to any large enough region for which the data of daily new cases are available. We have applied the SUTRA model to thirty-two countries, covering more than half of the world's population. Our conclusions are: (i) The model is able to capture the past trajectories very well. Moreover, the parameter values, which we can estimate robustly, help quantify the impact of changes in the pandemic characteristics. (ii) Unless the pandemic characteristics change significantly, the model has good predictive capability. (iii) Natural immunity provides significantly better protection against infection than the currently available vaccines.


Introduction
The COVID-19 pandemic caused by the SARS-CoV-2 virus has by now led to more than 600 million reported cases and more than six million deaths worldwide, as of October 1, 2022 [1]. By way of comparison, the infuenza epidemic of 1957 led to 20,000 deaths in the UK and 80,000 deaths in the USA, while the 1968 influenza pandemic led to 30,000 deaths in the UK and 100,000 deaths in the USA [2]. In contrast, the COVID-19 pandemic has already led to more than one million deaths in the USA and more than 170,000 deaths in the UK [2]. Therefore the COVID-19 pandemic is the most deadly since the Spanish Flu pandemic which started in 1918. In the USA, 675,000 people, or 0.64% of the population, died in that pandemic [3], compared to 0.33% of the population in the current pandemic. In order to cope with a health crisis of this magnitude, governments everywhere require accurate projections of the progress of the pandemic, both in space and over time, and at various levels of granularity. In addition, decision-makers also require assessments of the relative effectiveness of different non-pharmaceutical interventions (NPIs) such as lockdowns.
Over the past century or so, various mathematical models have been developed to predict the trajectory of a pandemic. These can be classified in three broad categories: statistical models, state-space models, and empirical models [4]. The most popular among these are compartment models, a type of state-space models. These models divide population into disjoint compartments representing different stages of infection, and are based on the premise that the disease spreads when an infected person comes into contact with a susceptible person.
In the initial SIR model [5], the population was divided into three compartments: S (Susceptible), I (Infected), and R (Removed / Recovered). Subsequently an intermediate compartment of E (Exposed) was introduced between S and I [6]. In SEIR model, interactions between S and E do not lead to fresh infections. For pandemics like COVID-19, that have significant number of asymptomatic patients, instead of E, compartment A (Asymptomatic) with interactions between S and A also leading to new infections, is more suited [7]. For COVID-19 pandemic, a number of models have been introduced to capture its trajectory. Two significant features of this pandemic are the presence of a large number of undetected cases and time-varying parameter values due to the emergence of various mutants, lockdowns etc. For any model to capture the trajectory well, it must take into account these two features. Some of the proposed models have large number of compartments in order to make the model biologically more realistic (for example [8]). These models have a large number of parameters, and estimating their values reliably is not possible from reported data due to the well-known phenomenon of the bias-variance tradeoff in statistics. This is even more true when the parameter values change with time. Besides the ones mentioned above, there are other models work with four compartments (and consequently a small number of parameters). In these models, parameter estimation is easier (for example [9,10,11]). However, even these models have to make certain assumptions about how the parameter values change. For example, [9] assumes the parameter values are derived from other considerations and do not change with time, while both [10] and [11] assume that parameters change following a specific equation.

Our Contributions
Against this background, in this paper we propose a four compartment model called SUTRA (Susceptible, Undetected, Tested positive, and Removed, Approach). 1 The components are same as in [10], however, there are differences in the dynamics (with better epidemiological justification, see 3). This has following consequences.

Fundamental Equation
The model admits a time invariant relationship (called fundamental equation) between detected new cases (N T ), detected active cases (T ), and total detected cases (C T ): whereβ andρ are parameters of the model and P 0 the population of region under study (see section 3).

Estimation ofβ andρ
We can efficiently estimate values of the two parametersβ andρ for the entire duration of the pandemic using the fundamental equation and simple linear regression (see section 5). These parameters have the following interpretation: • Parameterβ ≈ β, a standard parameter denoting contact rate (also called transmission rate).
• We may interpretρ ≈ where is detection ratio, the ratio of detected new cases to total number of new cases (as done in an earlier version of our model [12]), however, it is not satisfactory since in every region, value ofρ is observed to increase by a very large factor (1000 or more) in the initial couple of months of pandemic before becoming much less volatile (see section 9). A better interpretation isρ ≈ ρ where ρP 0 is effective population under the pandemic influence (see section 3). The large initial increase then makes sense since the effective population under pandemic at the beginning is a very small fraction of P 0 and increases very rapidly.

Phases of Pandemic
The value of parameterβ reduces when restrictive measures like lockdowns are imposed, and increases when these measures are lifted or a more infectious mutant arrives. Change in the value of parameterρ happens for many reasons. For example, when pandemic spreads to newer regions, or is completely eliminated from a region, or a part of susceptible population gets vaccine-induced immunity, or a part of population with immunity (acquired through vaccination or prior infection) loses it (see section 6). When the value ofβ orρ changes significantly, the trajectory of the pandemic changes. The model captures it as a phase change and recomputes the new values. As explained in section 6, the model can detect when the values ofβ orρ are changing, and when do they stabilize.

Future Projection
With the knowledge ofρ andβ, model can efficiently compute trajectory of the pandemic for the entire duration. The computed trajectory is a good estimate for future also as long as parameters do not change significantly (see section 7).

Estimation of ρ and
To understand the impact of pandemic better, it is desirable to estimate values of ρ and separately instead of their product. We show (Theorem 1) that given values ofβ andρ along with number of total and active infections on starting date of the simulation, there are only finitely many possible values for ρ and . Further, there is a unique canonical value for both at each time instant.
Since infections count at the start of simulation is not known, one requires a more realistic condition to be able to compute ρ and values. Towards this, we show that the condition can be replaced by knowledge of value of or ρ at any one time instant (section 8). Former requirement is met with a good serosurvey of the region at any point in time. Latter requirement is achievable if we can identify the time when pandemic has spread all over the region making ρ close to 1. One also requires that there is little vaccine-induced immunity at the time, since vaccine immunity reduces ρ (see Lemma 2).
For COVID-19, as the Omicron mutant arrived nearly eighteen months after pandemic started and is supposed to have bypassed vaccine-immunity nearly completely (as we also show in section 10), we can assume ρ ≈ 1 sometime after Omicron reached a region that did not implement strict control measures at the time. For such regions (that cover almost the entire world barring exceptions like China), we can estimate values of ρ and .

Analysis of the Past
The computed parameter values provide a quantification of impact of various events during the course of the pandemic. This includes impact of lockdowns and other restriction measures and arrival of new mutants. Section 9 does this analysis for four countries.

Analysis of Immunity Loss
The Omicron mutant caused widespread loss of immunity. Applying the model on thirty-two countries covering all continents and more than half the world's population, we deduce that loss of vaccine-immunity was significantly more than natural immunity conferred by prior exposure to any variant (see section 10).
This, coupled with the fact that vaccines continue to protect against severe infection, strongly suggests that the best strategy to manage the pandemic is to allow it to spread after vaccinating the population.

Model Formulation
Perhaps the earliest paper to propose a pandemic model incorporating asymptomatic patients is [7]. In this paper, the population is divided into four compartments: S, A (for Asymptomatic), I (for Infected) and R. Interactions between members of S and A, as well as between members of S and I, can lead to fresh infections. In that paper, it is assumed that almost all persons in A escape detection, while almost all persons in I are detected by the health authorities. While the SAIR model of [7] is a good starting point for modeling diseases with asymptomatic patients, and has been used in a few models for COVID-19 ([9] for example), it is not a good fit for COVID-19 for the following reasons: (i) due to contact tracing, some fraction of A does get detected and is often of similar order as detected symptomatic ones, (ii) many symptomatic cases are not detected. Therefore, the size of I cannot be estimated well.
In the present paper we propose a different grouping, namely: S = Susceptible Population, U = Undetected cases in the population, T = Tested Positive, either asymptomatic or symptomatic, and R = Removed, either through recovery or death. This leads to the SUTRA model, where the last A in SUTRA stands for "approach." Same division is used for models in [10,11]. As is standard, we use symbols S, U , T , R to also represent (time varying) fractional size of the four compartments.
The category R of removed can be further subdivided into R U denoting those who are removed from U , and R T denoting those who are removed from T . As in the conventional SAIR model [7], interactions between members of S on the and members of U or T , can lead to the person in S getting infected with a certain likelihood.  Figure 1. Typically, to handle undetected cases, models assume that the size of T is fraction of size of U +T (see for example, [10,11]). This is essentially equivalent to the assumption that detected new cases are fraction of new infections, as assumed in our model. Epidemiologically, all the new cases (but for rare exceptions) will remain undetected for a few days (until the symptoms appear). Therefore, one needs to justify the choice of βSU for detected new cases. We argue as follows: • Recently infected persons have higher chances of getting detected for two reasons. For symptomatic cases, the symptoms appears within a few days. For asymptomatic cases, they are detected through contact tracing which mostly starts with a symptomatic case and the asymptomatic cases detected would all be infected after the initiating symptomatic case.
• Number of new cases do not change dramatically over a few days and so number of detected cases over past few days can be taken to be proportional to βSU , number of most recent cases.
A few additional reasonable assumptions have been made to simplify parameter estimation. Specifically, • It is assumed that the removal rate for both compartments T and U is the same. This can be justified because, due to contact tracing, a significant fraction of patients in T are asymptomatic, and those people recover at the same rate as the asymptomatic people in U .
Even for the small fraction in T who develop complications and pass away, the time duration is very close to that of those who recover.
• There is no interaction shown between the T and S compartments. In most countries, those who test positive (whether symptomatic or not) are either kept in institutional quarantine, or told to self-quarantine. In reality, there might still be a small amount of contact between T and S. However, neglecting this does not significantly change the dynamics of the model, and greatly simplifies the parameter estimation.
With these considerations, the governing equations for the SUTRA model are: Since these quantities denote the fraction of the population within each compartment, we have There are three parameters in above equations, namely β, γ, and . The interpretation of these parameters is as follows: • β = The expected number of susceptible persons infected by an infected person in one day; it is called the contact rate or transmission rate.
• γ = Removal rate, the rate at which infected people are removed including both recoveries and deaths.
• = Rate at which infected patients in U move over to T . As shown later, it also equals the ratio T /(U + T ) most of the time, and is thus called the detection rate.
Later, we introduce two more parameters ρ and c, and derive expressions forβ andρ in terms of β, , ρ, and c.

Analyzing Model Equations
Defining M = U + T , R = R U + R T , we get from equations (2) and (3) thaṫ resulting in for an appropriate constant of integration c. Adding equations (2) giveṡ resulting in for some constant d. Since e −γt is a decaying exponential, it follows that, except for an initial transient period, the relationship M = 1 I holds. This in turn implies that U = M − T = 1− T . How long is the transient period? Observe that the constant d equals M (0) − 1 T (0) which is close to zero since fraction of infected cases at the start of pandemic is very small. Therefore the transient period will not last more than a few days. As we will see later, such transient periods will recur at various stages of pandemic and all of them remain small.
Define N T =Ṫ +Ṙ T = βSU , the fraction of population detected to be positive at time t, and C T = T + R T , the fraction of population detected to be infected up to time t. The above simplifications allow us to rewrite equation (2) as: Rearrange (8) as

Discretization of the Model Relationships
The progression of a pandemic is typically reported via two daily statistics: The number of people who test positive, and the number of people who are removed (including both recoveries and deaths). The second statistics has a problem though: there is no agreement on when to classify an infected person as removed. Some do it when RTPCR test is negative, some do it when symptoms are gone for a certain period, and some others do it after a fixed period of time. For the purpose of modeling, this classification needs to be done at the time when an infected person is no longer capable of infecting others. This is hard to decide, and so is almost never done. Further, some countries do not report second statistics at all (UK for example). In such a situation, we cannot rely on reported data, and instead compute R T by fixing γ to an appropriate value as discussed in section 4. Let T (t) denote the number of active detected cases on day t, R T (t) denote the number of detected cases that are removed on or before day t, and N T (t) denote the number of cases detected on day t. Note that all three are integers, and t is also a discrete counter. In contrast, in the SUTRA model, T , R T and N T are fractions in [0, 1], while t is a continuum. Therefore, where P is the effective population that is potentially affected by the pandemic. Now we introduce the parameter measuring the spread of the pandemic. Define number ρ, called the reach, which equals P/P 0 , where P is the effective population and P 0 is the total population of the group under study, e.g., the entire country, or an individual state, or a district (this parameter is also introduced and studied in [11]). The reach parameter ρ is usually nondecreasing, starts at 0, and increases towards 1 over time (situations where it decreases are discussed later). While the underlying population P 0 is known, the reach ρ is not known and must be inferred from the data.
Substituting P = ρP 0 , and integrating equation (9) over a day gives a relationship that involves only measurable and computable quantities T , C T = T + R T , and N T , and the parameters of the model, namely Note that N T is shifted forward by one day since new infections reported on day t+1 are determined by active infections and susceptible population on day t. Eq. (10) is the fundamental equation governing the pandemic. It establishes a linear relationship between N T , T , and C T T , that can be computed using the fundamental equation and the first three equations below (after fixing γ).
In addition to the fundamental equation, we will need discrete forms of other equations of the model to compute all quantities. We group them in two-first the quantities that can be computed from N T : The second group is of equations that involve numbers that cannot be computed from reported data: It is easy to see that all the quantities can be computed using above equations in addition to the fundamental equation once the parameter values used in the equations are available.

Fixing γ
As discussed in the previous section, reported removal data does not provide a good estimate for γ.
In [13], median duration of infection for asymptomatic cases was estimated in the range [6.5, 9.5] and mean duration for symptomatic cases in the range [10.9, 15.8] days with a caveat that the duration reduces when children are included. In [14], infection duration for symptomatic cases was observed to be less than 10 days. Since our groups U and T consist of a mix of asymptomatic and symptomatic cases, and it is likely that an infected person stops infecting others before becoming RTPCR negative, we take the mean duration of infection for both groups to be 10 days, implying γ = 0.1. All our simulations are done using the above value of γ and show a good fit with the actual trajectories.

Estimation ofβ andρ
One of the distinctive features of our approach is a methodology for estimating the values of all the parameters in the pandemic model from reported raw data on the number of daily new cases. The model has five parameters γ, β, , ρ, and c. At a first glance, these appear all independent, however, we show in section 8 that last four are essentially determined byβ andρ. In this section, we show how to estimateβ andρ from reported data using the fundamental equation. Let N T (t) be the reported new infections on day t. Note that N T (t) may not be the same as detected new infections on day t since there may be delays in reporting detected cases. Moreover, weekends often see fewer tests being done, causing unexpected variations in N T . To remove latter, we average N T (t) over a week, and let Let C T (t) = t s=0 N T (s), the total number of reported cases until day t, and T (t) be the number of reported active cases on day t computed inductively using equation Then the following linear regression problem is solved: The quality of the fit parameter, usually denoted by R 2 , is computed as follows: with the optimal parameter choices. The closer R 2 is to one, the better is the quality of the fit. At times, when there are relatively few data points (t 1 − t 0 is small), or the data has significant errors, above linear regression method fails to work (e.g., estimated parameter value becomes negative). In such situations we use a different method for estimation that is more tolerant to errors as described below. Let Find values ofβ > 0 andρ > 0 that maximize the product R 2 = R 2 β · R 2 ρ . This choice ensures that bothβ andρ play almost equally significant roles in minimizing the error. Further, the desired maximum of R 2 β R 2 ρ is guaranteed to exist: The only situation when the above method will not yield the desired maxima of R 2 is when u is dependent on either v or w. Former implies that T is proportional to N T over the time period, or equivalently, S does not change over the period. This implies N = 0 = N T = T for the period. Similarly, latter implies that T is proportional to C T T for the duration, or equivalently, C T does not change over the period. This also implies that N T = 0 = N . Either case occurs when the pandemic has effectively ended and there are no new cases for an extended period.
The uncertainty in the parameter estimation is computed using the standard mean-square error formula for linear regression. We use it to compute 95% confidence interval ranges forβ andρ values.

Phases of the Pandemic
The parameters ρ, β and are not constant, and vary over time. This causes changes inβ andρ as well. The contact rate β changes for following reasons: • Emergence of new and more infectious variants of the virus, which would spread faster than its predecessor. It takes time for the new variant to overtake whatever existed previously, which is why this factor would cause β to increase over a period.
• Non-compliance with COVID guidelines. The β parameter measures the likelihood of infection when an infected person (from either U or T ) meets a susceptible person from S. Thus β increases if people do not wear masks, or fail to maintain social distancing, and the like.
• The parameter can also decrease suddenly, with almost a step change, due to non-pharmaceutical interventions such as lockdowns.
The reach ρ changes for following reasons: • Spread of the pandemic to parts of the region that were previously untouched by it causes ρ to increase. The parts may even be physically co-located with parts already touched by the pandemic comprising of those people who had completely isolated themselves.
• Elimination of the pandemic from parts of the region that were under its influence causes ρ to decrease by the fraction of still susceptible population of the parts.
• Vaccination of susceptible people causes ρ to decrease, as these people moving out of susceptible compartment can be viewed as effective population under the pandemic reducing. Similarly, loss of immunity among immune population causes ρ to increase as this can be viewed as effective population under the pandemic increasing. This is formalized by the following lemma.
Lemma 2. Suppose ρ gain is the fraction of susceptible population that became immune via vaccination, and ρ loss is the fraction of immune population that lost immunity over a specified period of time. Then the new trajectory of the pandemic is obtained by multiplying both β and ρ (equivalently bothβ andρ) by 1 + Finally, the detection rate may increase due to more comprehensive testing, and may decrease due to reduction in testing.
The changes in parameter values occur either as a slow drift over an extended period of time, or as sudden rise and fall. We divide the entire timeline of the pandemic into phases, such that within each phase, the parameters are (nearly) constant. A phase change occurs when one or more parameter values change significantly. It could be due to a quick change for reasons listed above, or accumulated slow change over an extended period. By convention, we include the duration of change in a parameter as part of new phase and call it drift period of the phase. The remaining duration of a phase is called stable period of the phase.
When the value of changes, then the relationship T = M breaks down. The following lemma shows that T converges to M as soon as stabilizes to its new value. 4. Repeat the same with every subsequent phase.
We demonstrate the above methodology for one phase (phase #9) in India: when the deltavariant started spreading rapidly in the country during April 2021. In Appendix A, we have plotted During the drift period of a phase, when the parameter values are changing, the points continuously drift away from a line passing through the origin ( Figures 13,14,15) indicating that equation 10 is not satisfied. When the phase stabilizes, the points corresponding to the period line up nicely (Figures 16,17,18,19,20) indicating that the equation 10 is now satisfied. The plots also show that values ofβ andρ are changing quickly during the drift period, and do not change much during stable period. This leads to easy identification of phases and stable period within.

Parameter values during drift period
We have so far seen how to estimate values ofβ andρ during stable period of every phase. However, in order to simulate the course of the pandemic, it is necessary to have the values of the parameters during the drift period as well.
Suppose d is the number of days in drift period, and b 0 and b 1 are the computed values of a parameter in the previous and the current phases. Then its value will move from b 0 to b 1 during the drift period. A natural way of fixing its value during the period is to use either arithmetic or geometric progression. That is, on ith day in the drift period the value is set Among these, geometric progression captures the way parameters change better: • When a new, more infectious, mutant spreads in a population, its infections grow exponentially initially. This corresponds to a multiplicative increase in β.
• Similarly, a new virus spreads in a region exponentially at the beginning. This corresponds to a multiplicative increase in ρ.
• A lockdown typically restricts movement sharply causing a multiplicative decrease in β.
• A change in testing strategy typically gets implement fast in a region, causing a multiplicative change in .
For these reasons, we assume that changes in parameters β(1 − ) (this is the effective contact rate due to quarantining of detected cases), ρ, and are multiplicative. Further, changes in parameter c are additive as it is constant of integration ensuring continuity between two phases. Therefore, we may assume that changes in 1 − c are multiplicative. This leads to the conclusion that changes inβ andρ are also multiplicative.
Having defined how the parameters change during drift periods, we assume that the equations (10), (11), and (12) hold on all days. When in drift period, even the parameter values in the equations change daily as defined above. Subsequent sections show that our model with these assumptions is able to capture the trajectory of the pandemic very well.

Future Projections
Once the quantitiesβ,ρ are estimated as above for current phase, equations (10) and (11) can be used to compute values of N T , C T , and T for the entire phase duration. If the model captures the dynamics well, the predictions for daily new cases N T should match closely with averaged reported numbers N T after the phase enters stable period as long as parameters do not change significantly. Indeed, this is confirmed by our simulations of trajectories in multiple countries. For example, for the phase #9 of India discussed in the previous section, predicted trajectory changed rapidly when the phase was in drift period, and stabilized when it transitioned to stable period (see Figure 2). This property allows one to accurately predict the future course of the pandemic once the present phase stabilizes. We used this to make several successful predictions in the past. Some notable ones were predicting the timing and height of the peak of second wave of India ten days in advance [15], predicting timing of the peak of third wave in India as well as many states of the country [16], predicting timing and height of the peak of Delta-wave in UK ten days in advance [17], and predicting timing and height of the peak of Delta-wave in US more than a month in advance [18]. The predictions for India and its states were useful to the policy-makers in planning the required capacity for providing health care, and scheduling nonpharmaceutical interventions such as school reopenings.

Estimation of ρ and
After fixing values of γ, the model has four parameters left: β, ρ, and c. We have seen how to estimate values of composite parametersβ andρ at all times, which allows us to compute the trajectory of daily new detected cases for a region. We can obtain more information about the pandemic if values of ρ and can be estimated separately. For example, will enable us to estimate trajectory of total daily new cases, including undetected ones. We provide more applications in the next two sections.
Without any additional information, besides the daily new detected infections time series, it is not possible to estimate value of : In this section, we show that with just one additional data point-C(0) and M(0), the total number of cases up to time t = 0 and total active cases at t = 0--the number of possible trajectories for N (t), consistent with the given data, becomes finite: Theorem 1. Given detected new cases trajectory, N T (t), 0 ≤ t ≤ t F and C(0), there exist only finitely many trajectories for N (t) consistent with N T . Further, a good estimate for all the trajectories can be obtained efficiently.
As is shown in the proof of above theorem (see Appendix C), the trajectory for N (t) is unique for the first phase, but there may be multiple ones for subsequent phases, identified by a unique value for the pair ( , c) for each. Using the observation that the value of from one phase to next will not change significantly, we can identify a unique canonical trajectory for total new cases: for each phase, given the possible values of that give rise to consistent trajectories for N (t), choose the canonical value of to be the one closest to the canonical value of of previous phase (for first phase, there is anyway a unique value of ). The corresponding trajectory for N (t) is called canonical trajectory.
As the theorem also states, the canonical value of and corresponding value of c can be efficiently estimated. This, in turn, provides values of β =β (1− )(1−c) and ρ =ρ (1−c) for the phase.
In this way, we get the values of all parameters at all times.

Calibrating the Model
Above shows how to estimate parameter values for all phases, provided we know the values of C(0). This is equivalent to finding out the values of parameters for the first phase, say 1 , since C(0) = 1 1 C T (0) (since c 1 = 0 as shown in the proof of theorem 1). While this is good in theory, we do not know 1 in practice. Moreover, time t = 0 when the data becomes available for the first time is unlikely to be the time when the pandemic begins, and hence we may not have c 1 = 0. However, c 1 = 0 will still be a good estimate since R(0) is likely to be very small. We estimate value of 1 from other information available about the pandemic. This is called calibrating the model. We can calibrate the model in two ways: • A sero-survey at time t 0 provides a good estimate of C(t 0 − δ), where δ equals the time taken for antibodies to develop. Once we accurately estimate 1 , the model can compute C(t) at all times t. We choose a suitable value of 1 ensuring that model computation matches with the sero-survey result at time t − δ.
• When the pandemic has been active long enough in a region without major, long-term restrictions, we may assume that it has reached all sections of society, making ρ close to 1. Again, we can choose 1 that ensures that the reach of the pandemic is close to 1 at suitable time.
While using the above two methods for calibrating the model, following points need to be kept in mind: Using serosurveys. Many serosurveys suffer from significant sampling biases. For example, if a survey is done using residual sera from a period of high infection numbers, it is likely to significantly overestimate the seroprevalence because a large fraction of uninfected persons would not venture to give blood sample in such a period. In order to minimize sampling biases, therefore, one should use serosurveys done during a period of low infection numbers. Even then, some uninfected people may not participate making the estimates higher than actual. To further reduce bias, one should ideally be able to use multiple serosurveys as well as use the fact that reach is close to 1 by a given time.
Using reach. As observed earlier (Lemma 2), parameter ρ is impacted by several factors, including gain and loss of immunity. Therefore, ρ may not be close to 1 even when the pandemic has spread over entire population. To capture this, we define ρ actual to denote the actual reach of pandemic, so ρ = ρ actual · (1 + ρ loss −ρ gain ρ actual ) = ρ actual + ρ loss − ρ gain . To calibrate the model using ρ actual ≈ 1 at certain time, we need an estimate of ρ loss − ρ gain too, which introduces more errors in calibration.
There are regions where neither an accurate sero-survey is available, and it is evident that reach is nowhere close to 1. For such regions, calibration cannot be done with any confidence, and so estimation of all parameter values is not possible.

Analysis of the Past
The parameter table of a country enables us to quantify the impacts of various events like the arrival of a new mutant, or a lockdown. Moreover, through the reach parameter, we can also explain the somewhat mysterious phenomenon of multiple peaks occurring in rapid succession that was observed in many countries.
Below, we discuss in detail the progression of the pandemic in four countries. The time series data for India was sourced from [19], and for rest of the countries from [1].

India
Model computed trajectory for detected cases shows an excellent match with the reported trajectory: For estimating all parameters, the calibration was done using sero survey done in December 2020 [20], a period of low infection. Estimates of seropositivity computed by the model were matched with two other serosurveys [21,22], and very good agreement was found. Further, our model showed that the reach was close to maximum by December 2021, a very likely scenario. Interestingly, the detection rate stayed almost unchanged at 1/32 throughout the course of the pandemic. Comparing the timeline of the pandemic [23] with parameter table below, we observe the following.
First wave (March to October 2020): The strict lockdown imposed at the end of March 2020 brought down the contact rate β by a factor of two. The reach was very small until May (≈ 0.03) but increased to 0.36 between the end of June and the end of August. This was caused by reverse migration of workers and a partial lifting of lockdown that happened during this period.
Second wave (February to July 2021): The arrival of the Delta variant caused the value of β to rise to 0.38 in February 2021. As the variant began to spread in different parts of the country, most states imposed restrictions, which reduced the nationwide β to 0.28 by April.
In the same month, ρ increased sharply to 0.83.
Note that while increase in ρ was clearly due to Delta variant, the increase happened more than a month after increase in β. We have observed this delayed increase phenomenon in ρ repeatedly.
The removal of all restrictions by August caused β to increase to 0.58. This suggests that the Delta variant was more infectious by a factor of ≈ 2 compared to original variant.  A ripple (April 2022 to September 2022): The value of β increased to around 3 by June. This, coupled with an increase in ρ from 1.02 to 1.08 (indicating around 6% population losing natural immunity) caused a ripple that peaked in July.
At present, around 98% of population is estimated to have natural immunity.

UK
Model computed trajectory for detected cases shows a good match with the reported trajectory:

Figure 4: Predicted and Actual Trajectories for UK
For estimating all parameters, the calibration was done using two of the three serosurveys reported in [24]. During the period of serosurveys, the infection numbers were going up and down, which made calibration a little tricky. We used the numbers from the last two surveys as well as the observation that reach has remained stationary from November 2021 (suggesting that ρ actual has been around 1 since then) for calibration. The detection rate started at 1/9.3 and over time increased to almost one in three cases. Comparing the timeline of the pandemic [25] with parameter table below, we observe the following.
First wave (March to July 2020): The strict lockdown imposed in March 2020 brought down the contact rate β from 0.26 to 0.15 in mid-April. However, almost simultaneously, ρ increased three-fold causing another peak. By July, β was back up to 0.24 after removal of restrictions.
Second wave (September 2020 to January 2021): This wave was primarily caused by increase in value of ρ from 0.08 to 0.6. This increase in ρ was a natural consequence of very small effective population until August (less than one percent) and easing of lockdown from July (as noted above, increase in ρ happens with a lag). As the numbers started increasing, fresh restrictions were put in place bringing β down by 20% by September-end. This caused the cases to peak by October-end (by that time ρ increased to 0.3). As the lockdown was eased before Christmas, both β and ρ started increasing causing and second bigger peak in January 2021.
The rise in case number caused another lockdown, but this time β did not decrease. Note that a new variant, called Alpha, started spreading in UK rapidly in December 2020. It was believed to be significantly more infectious than earlier one. This appears to be the reason why the value of β did not decrease in January, and went up slightly instead.
Third wave (February 2021 to October 2021): Lockdown was eased during February-March which resulted in a significant rise in β to 0.68 by second half of March. This jump, however, caused only a slight change in trajectory because reach stayed around 0.6 and more than 85% of population within reach had natural immunity by then. Numbers started rising from mid-June due to increase in ρ (again a delayed increase). There were three peaks in quick succession: The first caused by increase in ρ to 0.76 in July, the second caused by further increase in ρ to 0.92 in August (when β came down to 0.37 during this period, likely caused by precautions taken by people due to high numbers), and the third caused by increase in β to 0.5 in addition to a slight increase in ρ. This increase in β was likely due to the Delta variant now active in the country. At present, around 92% of population is estimated to have natural immunity.

US
Model computed trajectory for detected cases shows an excellent match with the reported trajectory:

Figure 5: Predicted and Actual Trajectories for USA
For estimating all parameters, the calibration was done using the serosurvey [26]. The samples were taken from life insurance applications. The calibration was further supported by the fact that ρ has not changed since December, suggesting that ρ actual was close to 1 at the time. The detection rate has slowly decreased from 1/3.5 to 1/4 during the course of the pandemic. Comparing the timeline of the pandemic [27] with parameter table below, we observe the following.
First wave (March to August 2020): Restrictions imposed in April 2020 brought down the contact rate β from 0.31 to 0.18 by mid-May. However, almost simultaneously, ρ increased to 0.04 causing a flat trajectory. In June, most restrictions were lifted. This increased ρ further to 0.11 causing a peak in July-end. The value of β, however, did not increase. This could be due to precautions taken by a large number of people.
Second wave (September 2020 to February 2021): By October, β went up to 0.24 and stayed around this value until the end of the wave. The value of ρ increased in three steps: to 0.26 during September-October period, to 0.33 in December, and to 0.39 in January. This causes three successive peaks in November, December, and January.
Third wave (March 2021 to November 2021): There were two peaks separated by more than four months in this period. The Delta variant appeared to have arrived in March causing β to increase to 0.45. However, it caused only a small peak since ρ stayed around 0.5 until July, and more than 75% of population under reach had natural immunity. The reach started increasing in August to eventually become 0.7 by mid-September causing another peak (yet another case of delayed increase in ρ). December causing β to increase to 0.6, but the numbers did not increase much by Decemberend, since ρ did not change by much. Then the reach increased substantially to 1.08 in a short time leading to a very sharp and high peak. By February, the wave subsided, and even though β jumped to 1.87 in March, it did not cause cases to increase as more than 90% of population was immune by then.
Fifth wave (April 2022 to September 2022): The primary cause of this wave appears to be a loss of natural immunity. By May, ρ increased to 1.2 and is close to 1.3 at present. This implies that more than 20% of population has lost natural immunity in past six months.
Immunity loss at such a scale has not been observed in the other three countries discussed here. Reasons for this are not clear.
At present, around 85% of population is estimated to have natural immunity.

South Africa
Model computed trajectory for detected cases shows an excellent match with the reported trajectory:

Figure 6: Predicted and Actual Trajectories for South Africa
For estimating all parameters, the calibration was done using the serosurvey [28]. The calibration was further supported by the fact that ρ has not changed since November suggesting that ρ actual has been close to 1 since then. The detection rate has remained almost unchanged at 1/17 during the course of the pandemic. Comparing the timeline of the pandemic [29] with parameter table below, we observe the following.
First wave (April to August 2020): Restrictions imposed in March and April 2020 brought down the contact rate β to around 0.2. A significant increase in ρ to 0.24 by June-end caused the first wave that peaked in July-end.
Second wave (September 2020 to February 2021): Restrictions were lowered in September causing increase in β value to 0.25. Reach also continued to increase slowly to 0.5. This caused only a slow rise since more than 40% of population under the reach was already immune. The Beta variant arrived in December causing an immediate jump in β value to 0.46. This caused the second peak in January. Restrictions were reimposed in December to control the rise in the numbers due to Beta variant. The fact that β still increased substantially shows that infectiousness of this variant was quite high.
Third wave (March 2021 to October 2021): With removal of restrictions measures by March, value of β further increased to 0.65. However, since increase in ρ is typically delayed and more than 90% of population within reach was already immune, the rise in β did not cause increase in numbers. The numbers started rising when ρ started increasing in April to become 0.75 by mid-May. Restrictions were partly brought back causing β to come down to around 0.4. The value of ρ further went up to 0.94 causing a peak in July. As the numbers started coming down from the peak in August, an unusual phenomenon occurred. Cases started increasing once again, there was a short peak in second half of August, and then the numbers came down once again but with a slightly less steep slope than before. Our model shows that this happened due to a sharp increase in value of β to nearly 1 from 0.43. Note that restrictions were being increased during June-July and were relaxed only from September, so the increase in β was not due to relaxations. Was this caused by Omicron variant that was detected later in South Africa? The sharp increase in β which then stayed high certainly suggests so.
Fourth wave (November 2021 to March 2022): In November, with more and more relaxations, Omicron caused β to further increase to 1.58 and ρ to 1.03 resulting in a high peak. No restrictions were imposed this time, and so the numbers rose and fell sharply.
Fifth wave (March 2022 to June 2022): A further increase in β to around 3 by April resulted in another peak in mid-May. This peak was, however, a small one since reach was stationary around 1.05 and more than 95% of population had natural immunity.
At present, around 98% of population is estimated to have natural immunity.

Analysis of the Immunity Loss
After the South African authorities announced the emergence of a new variant of concern (VOC), later named Omicron, the epidemiology community started analysing the ability of the Omicron variant to bypass immunity provided by vaccination, or prior exposure, or both. Our objective in this section is to provide a quantitative analysis using the SUTRA model. But before that, we give a brief summary of the vast literature based on laboratory (as opposed to population-level) studies. Everywhere in the world where it was discovered, the Omicron VOC soon replaced all other variants and was responsible for a massive increase in cases. This was due to high transmissibility conferred by the mutation, ensuring a tight binding to the ACE 2 receptor facilitating immune escape [30]. The immune escape phenomenon was reported by many groups studying the neutralization activity of sera from both infected and vaccinated individuals; see [31,32,33,34,35]. The immunity conferred by complete vaccination decreased from 80% for the Delta variant to about 30% for the Omicron variant. People infected with the Delta were better off than those infected with the initial Beta variant. There was a complete loss of neutralizing antibodies in over 50% of the vaccinated individuals and the decrease in titres varied from 43-122 fold between vaccines [36]. A booster Pfizer dose could generate an anti-Omicron neutralizing response, but titres were 6-23 fold lower than those for Delta variant. Sera from vaccinated individual of the Pfizer or Astra Zeneca vaccine barely inhibited the Omicron variant five months after complete vaccination [37]. In addition, Omicron was completely or partially resistant to neutralization by all monoclonal antibodies tested [30]. Overall, most studies confirmed that sera from convalescent as well as fully vaccinated individuals irrespective of the vaccine (BNT162b2, mRNA-1273, Ad26.COV2.5 or ChAdOx1-nCoV19, Sputnik V or BBIBP-CorV) contained very low to undetectable levels of nAbs against Omicron. A booster with a third dose of mRNA vaccine appeared to restore neutralizing activity but the duration over which this effect may last has not been confirmed. Double vaccination followed by Delta breakthrough infection, or prior infection followed by mRNA vaccine double vaccination, appear to generate increased protective levels of neutralizing antibodies [38]. Viral escape from neutralising antibodies can facilitate breakthrough infections in vaccinated and convalescent individuals; however, pre-existing cellular and innate immunity could protect from severe disease [38,39]. Mutations in Omicron can knock out or substantially reduce neutralization by most of the large panel of potent monoclonal antibodies and antibodies under commercial development. Studies also showed that neutralizing antibody titers against BA.2 were similar to those against the BA.1 variant. A third dose of the vaccine was needed for induction of consistent neutralizing antibody titers against either the BA.1 or BA.2.3,4 variants, suggesting a substantial degree of cross-reactive natural immunity [40].
The studies above indicate that vaccine immunity was lost substantially against Omicron, and natural immunity provided better protection. All the studies were done in laboratories or in a small section of population, and our analysis in this section complements them as it is based on population-wise data.

Vaccine Immunity before Omicron
We first analyze the gain in immunity due to vaccination before the arrival of Omicron. For this propose, we downloaded an extensive list of serosurveys, carried out in various countries and maintained by the site [41], eliminated surveys that were not done at national level, or had small sample sizes, or had high risk of bias. Nineteen countries remained after this pruning. These sero-surveys were used together with the SUTRA model to capture the pandemic trajectories and estimate parameter values in these countries. We identify two values for each country: 1. Value 1 ρ − 1 before arrival of Omicron. All the nineteen countries had restrictions removed well before Omicron and therefore, it is reasonable to expect that ρ actual was close to 1 by the arrival of Omicron in the country. As shown in Lemma 2, ρ = ρ actual + ρ loss − ρ gain ≈ 1+ρ loss −ρ gain at the time. In other words, ρ gain −ρ loss ≈ 1−ρ. Since calibration of the model provides only an approximate value of ρ, we use fractional gain in immunity ρ gain −ρ loss ρ ≈ 1 ρ −1 which is likely to be more robust.
2. Fraction of uninfected population in the country that has received at least one dose of vaccination at the onset of Omicron wave. To estimate this number, we assume that the two types of immunity, vaccine and natural, are independent random variables, implying that the fraction with hybrid immunity is the product of vaccine immunity and natural immunity fractions. With this assumption, and using vaccination data from [42], we can estimate the required value. Figure 7 plots the above two numbers. It shows a very strong correlation between the two numbers implying that vaccination provided excellent immunity before Omicron.

Loss of Vaccine Immunity after Omicron
We can measure immunity loss due to Omicron by comparing the value of ρ after Omicron arrives in a country with the value before its arrival. This change will be almost entirely due to immunity loss since ρ actual ≈ 1 before Omicron as discussed above.
We first compare it with vaccine-only immunity present in the population to get an estimate of how much of it was lost. Figure 8 plots these two numbers. Again, a very strong correlation is observed between the numbers. This, and the fact that the slope of best-fit line is close to 1, suggests that almost all of immunity loss was due to loss of vaccination immunity.
The above conclusion is further strengthened by the next plot where we compare immunity loss due to Omicron with the natural immunity present in the population before the arrival of mutation. Figure 9 plots these two numbers. It shows a very strong negative correlation implying that natural immunity provided excellent protection against Omicron.

Incorporating More Countries
To make our conclusions more broad-based, we include seventeen more countries based on following criteria: 1. All continents are represented well (five from Africa, two from North America, four from South America, thirteen from Asia, eleven from Europe, and one from Australia) 2. Populous countries are simulated (except China for which it is not possible to calibrate the model). More than half the world's population lives in these countries.
3. It is likely that ρ actual was close to maximum in these countries at the time of Omicron's arrival, allowing us to calibrate the model.
Adding these countries to the plots, we find little change in the correlations (see Figures 10,11,12), further strengthening the conclusions.
Taken together, these plots show that Omicron bypassed vaccination immunity almost completely, but natural immunity provided excellent protection. A clear conclusion is that countries that followed zero-COVID strategy -strictly control the spread and vaccinate entire populationsuffered maximum during the Omicron wave. Indeed, a perusal of Table 5 shows that in countries where the reach was very low before the arrival of the Omicron variant saw very large percentage increases in reach thereafter. It also suggests that the best strategy for managing the pandemic for a country is to allow the virus to freely spread after vaccinating entire population. Trying to control its spread even after vaccination will not build natural immunity in the population and there will always be a chance of fresh outbreaks. We can see it happening in China at present.

A India Phase 9 Plots
All the plots in this section have 1 P 0 C T T on x-axis and T − 1 β N T on y-axis.

B Pandemic Status in Countries at the Time of Omicron Arrival
In this section, we provide the data used for analysis in section 10. The table below lists, for 32 countries, the percentage of population that had vaccine immunity (given at least one shot of vaccine fifteen days before) when Omicron arrived in the country. It also lists percentage of population with natural immunity (estimated by the model), with hybrid immunity (assuming two types of immunity are independent), with vaccine-only immunity (difference vaccine immunity and hybrid immunity), and within reach of the pandemic at the time.

C Proofs
In this Appendix, we provide proofs of all lemmas and theorems stated in the paper.
Lemma 1 When u is independent of v as well as w, there is a maxima of R 2 with R 2 β , R 2 ρ ,β,ρ > 0.
Proof. Let x = 1 β and y = 1 ρ . Then we have: The denominator of equation (13) is always positive since u is independent of v as well as w. The numerator is a product of four linear terms in the unknowns x and y. Therefore the value of R 2 is positive inside the polygon defined by: and is zero on the boundaries. This guarantees that there exists at least one maxima inside the polygon.
Lemma 2 Suppose ρ gain is the fraction of susceptible population that became immune via vaccination, and ρ loss is the fraction of immune population that lost immunity over a specified period of time. Then the new trajectory of the pandemic is obtained by multiplying both β and ρ (equivalently bothβ andρ) by 1 + ρ loss −ρ gain ρ .
Proof. During the course of the pandemic, first few phases did not have any vaccination or immunity loss. Consider the first phase with either immunity loss or gain through vaccination or both. Suppose fraction of population that gains immunity through vaccination in this phase is ρ gain and the fraction of population that loses immunity is ρ loss . Consider a time instant t in the stable period of the phase when the changes in immunity have already taken place. Then the fraction of removed population would be R(t) + ρ gain P 0 − ρ loss P 0 where R(t) is the fraction of removed population if there was no change in immunity levels. Therefore, the fundamental equation (10) changes to: . Therefore, fundamental equation now holds with values of β and ρ multiplied by f . The other equations of the model (11 and 12) can easily be seen to hold with the same change in values of β and ρ.
Lemma 3 Suppose a new phase begins at time t 0 with a drift period of d days. Further, suppose that the value of parameter changes from 0 to 1 during the new phase. Then, Proof. The model parameters β(1 − ), , ρ and 1 − c change multiplicatively until they stabilize to new values. Suppose that in one day, β(1 − ) changes by a factor of f b , by a factor of f e , ρ by a factor of f r and 1 − c by a factor of f c . Then, composite parameterβ = β(1 − )(1 − c) will change by a factor of f b f c andρ = ρ(1 − c) will change by a factor of f e f r f c . We know that after the drift period,β andρ stabilize. Suppose that continues to change even after the drift period of d days and stabilizes after D days in the phase. We now consider following cases: •β changes during the drift period and β(1 − ) changes on day D − 1. In that case, change iñ β on day D − 1 equals f b f c which is not equal to 1 sinceβ changes during drift period. This is not possible.  (15), as in the proof of Theorem 1, can be carried out for day D (since all parameters stabilize by then), and substituting y = 1 fc = 1 in the equation (15) we get f e = x = 1 (note thatρ D =ρ D−1 ). This contradictions the assumption that changes during the phase.
•ρ changes during drift period and ρ changes on day D − 1. Then change inρ on day D − 1 equals f r f e f c which is not equal to 1 sinceρ changes during drift period. This is not possible.
•ρ changes during drift period and ρ does not change on day D − 1. Then change inρ on day D − 1 equals f e f c which must be equal to 1 sinceρ does not change after d days. Therefore, f e = 1/f c . Going back to equation (15) and substituting x = f e = 1/f c = y, we again get f e = 1 = f c , contradicting the assumption that changes during the phase.
Together, the cases above cover all possibilities and hence we conclude that D = d and therefore, Lemma 4 Given detected new cases trajectory, N T (t), 0 ≤ t ≤ t F , there exist infinitely many total new cases trajectories and corresponding values of consistent with N T .
Proof. Given N T (t), 0 ≤ t ≤ t F , we can compute phases of the trajectory and values ofβ andρ for all phases as shown in section 5, as well as C T (t), T (t) and R T (t) for the entire duration. Choose any value 0 in the range [0.9, 1.0]. Fix = 0 and c = 0. This allows us to compute the values of β and ρ for all phases (the value of ρ will be at most 1 0.9 times the value ofρ at any time).
Let N (t) = 1 N T (t), and M(t) = 1 T (t) for 0 ≤ t ≤ t F . Setting R(0) = 0, and using the equation (12) for R, we get that R(t) = γ t−1 s=0 M(s) = γ t−1 s=0 T (s) = 1 R T (t). Then, for all t, 0 ≤ t ≤ t F : It is straightforward to see that the remaining model equations (12) are also satisfied. As there are infinitely many values in the range [0.9, 1.0], the proof is complete.
Theorem 1 Given detected new cases trajectory, N T (t), 0 ≤ t ≤ t F and C(0), there exist only finitely many trajectories for N (t) consistent with N T . Further, a good estimate for all the trajectories can be obtained efficiently.
Proof. Proof is by induction on the number of phases. In the base case we have only one phase. For this phase, there is no drift period since there are no previous values of parameters. Therefore, parameter values stay the same throughout the phase duration. Let β 1 , ρ 1 , 1 and c 1 be the parameter values governing the actual trajectory for this phase. Therefore, M(t) = 1 1 T (t) and R(t) = 1 1 R T (t) + c 1 ρ 1 P 0 for the entire phase. Note that R(0) = 0 since at time t = 0, when the pandemic starts, there are no recoveries. Hence, c 1 = 0. Further, 1 = T (0)/M(0) = T (0)/C(0). From this, we can compute giving values of all parameters for first phase, using which the trajectory can be computed for the first phase uniquely. Suppose there are finitely trajectories up to phase i − 1. Fix any one trajectory with values of four parameters in phase i − 1 being β 0 , ρ 0 , 0 and c 0 . Let t 0 be the time when phase i starts. We have: Suppose phase i has a drift period of d days. In the model, the parameter values change multiplicatively during the drift period. Let j = 0 x j , and 1 − c j = (1 − c 0 )/y j for 1 ≤ j ≤ d, where x and y are unknown multipliers by which the two parameters change every day. The final value of the parameters will be d = 0 x d and 1 − c d = (1 − c 0 )/y d .
Let β j =β j (1− j )(1−c j ) and ρ j =ρ j j (1−c j ) for 1 ≤ j ≤ d. Then we can write: M(t 0 + j) R(t 0 + j) = g j 0 γ 1 · M(t 0 + j − 1) Therefore, both M(t 0 + j) and R(t 0 + j) are polynomials in x and y. It is straightforward to show that the degrees of M(t 0 + j) and R(t 0 + j) equal 2 j − j − 1 and 2 j−1 − j − 2 respectively. At the end of drift period, we have: For what values of unknown multipliers x and y are the relationships in equation (14) satisfied? To see this, we analyze the quantities N (t 0 + d), M(t 0 + d), and C(t 0 + d). From (12) we have: Therefore, Since where the second equality is from fundamental equation (10), we have: Therefore, .
Equation (15) 1]. Further, there is a simple and efficient way to list out good estimates for all the solutions: since −1 < c d < 1, one can step through possible values of c d in the range using small discrete steps, compute the value of y for the chosen value of c d , compute value of x using equation (15), and then check if P (y) is close to zero. This will list out good estimates of all feasible solutions.