1.1 Disease progression

SARS-CoV-2 infection manifests as an acute respiratory infection, progressing to respiratory failure in a small number of patients [54–56]. It can result in a range of clinical manifestations, from asymptomatic or mild infection to severe, requiring hospitalization [57]. Among patients who are symptomatic, the median incubation period is approximately 4 to 5 days [58]. About 97% have symptoms within 11 days after infection [59]. These can further be categorized as mild, severe and critical [60]. Patients who are hospitalized, the severe category, can progress to severe pneumonia and acute respiratory distress syndrome (ARDS). A fraction of these patients will require ventilation and an even smaller fraction may die [61]. Older patients experience greater clinical severity of COVID-19, which we account for via age-stratified parameters [61, 62]. Co-morbidities such as cardiovascular disease, diabetes and obesity are common underlying conditions associated with worse clinical outcomes and increased disease severity [62–64]. These can be incorporated through a composite risk score affecting branching rates between mild and severe disease states, although we do not do so here. Males may experience more severe disease than females, and genetic variations, including the ABO blood type, have been implicated in clinical outcomes for patients with COVID-19. Our model ignores these effects [65, 66]. Around 40–75% of infections may be asymptomatic, a fairly broad range. We note that numbers for India suggest a larger fraction of asymptomatic cases than reported elsewhere [67, 68].

1.2 Transmission dynamics

The compartmental structure we consider contains susceptible (S), exposed (E), asymptomatic infectious (I^a), pre-symptomatic infectious (I^p), mildly symptomatic infectious (I^m), severely symptomatic infectious (I^s), hospitalized (H), dead (D) and recovered (R) compartments. Transitions between these model compartments are shown in Fig 2. The dotted lines indicate the force of infection from the infected compartments on the susceptible population. Frequency-dependent transmission is assumed. As appropriate for a fast-spreading infection, we ignore the effects of demography [51].

The model equations, for unstructured compartments, are (1) The left-hand side of these equations denote first-order time derivatives. The population size is N. Infectious individuals in any associated compartment can infect the susceptible population regardless of symptoms and severity with a fixed transmission rate β. However, this quantity is modulated by the relative intensity of contacts between susceptible and infectious individuals whose effect is simply specified here through factors of ϵ. These can, in principle, vary with time and can also be chosen to vary between compartments.

Each infected individual enters an ‘exposed’ compartment (E), spending an average latent period 1/γ days before becoming infectious [69]. A fraction α of this infectious population remains asymptomatic (I^a) until recovery. Asymptomatic individuals are assumed to be infectious for an average period of 1/λ_a days [70, 71]. Those in the remaining fraction, of (1 − α), enter a pre-symptomatic state (I^p) where clinical symptoms are not exhibited for a relatively short average period of 1/λ_p [70]. Individuals in this pre-symptomatic compartment go on to developing mild or severe symptoms [72]. All rates applicable to the model are provided in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Age dependent branching ratios between compartments.

Our values for α_i and μ_i are from those consolidated by the COVASIM program, Refs. [76] and [77]. We use IFR’s from data from LMIC’s taken from fits in the paper of Ref. [78], where the following formula is derived: log₁₀(IFR) = −3.27 + 0.0524 * age. We modify only the δ_i for the 80+ age group to account for the leveling off of mortality in older age groups described by Ref. [24]. The IFR for each age group can be obtained as IFR_i = (1 − α_i)(1 − μ_i)δ_i. We use the estimated δ_i numbers as our initial IFRs, although when we allow them to vary as part of our minimization strategy to find the effective IFR’s we multiply all δ_i’s by a single smooth time-varying factor, optimizing this against data.

https://doi.org/10.1371/journal.pcbi.1010632.t001

A fraction μ of symptomatic cases is assumed to develop mild symptoms (I^m), while the remaining fraction (1 − μ) of cases are transferred to the severe class (I^s). Infectious cases with mild symptoms recover without hospitalization after 1/λ_m days. Severe cases require hospitalization after an average of 1/λ_s days. From the hospitalized population (H), we assume that a proportion (1 − δ) recovers successfully (R), after spending an average duration of hospitalization 1/ρ days. The remaining fraction, δ, of the hospitalised, will die. The numerical values of these parameters are specified in Table 1 [73, 74].

Of these parameters, the infectivity parameter β is particularly central. It determines the effective reproduction ratio as the epidemic proceeds. The force of infection arising from asymptomatic cases alone is assumed to be lower in comparison to that arising from the pre-symptomatic, mildly symptomatic and severely symptomatic cases.

1.3 Time-dependence of effective contacts to model non-pharmaceutical interventions

The effective infectivity in the INDSCI-SIM model is a product of three terms. The first is a bare infectivity parameter (β_ij ≡ β)and the second is the term involving the contacts, the C_ij’s. The third is the temporal modulation, the factors of ϵ. While it is only variation in this overall product that is of significance, it allows us to conceptualize interventions in a more targeted way: the sudden changes in infectivity imposed by a lockdown can be associated with abrupt changes in β whereas slow and secular improvements in masking etc. can be associated to the smoothly varying factors of ϵ. This product is affected by mask-wearing, hand-washing, voluntary self-isolation or self-quarantining and the maintenance of social-distancing. It is also influenced by global non-pharmaceutical interventions such as closures of schools, malls and cinemas, large-scale lock-downs, restrictions on movement and autonomous modifications of social behaviour.

We first use a “global” value of β, one that differs in the lockdown period and in the open period. We then model compliance with restrictions by modulating this with a time-dependent factor, using a hyperbolic tangent function of time, as described below. This function comes with an associated time-scale which is an output of our minimization procedure. A quantification of mask-wearing India-wide comes from studies incorporated into the IHME model for India that are standardized using survey results from the University of Maryland Social Data Science Center, from the Kaiser Family Foundation and the YouGov COVID-19 Behaviour Tracker survey [79]. This data motivates the use of such an interpolating function, since it shows an initial lag period, a sharp rise as public awareness increases and saturation at about 70% reflecting social acceptance of mask-wearing.

We also allow for a decay of the relative intensity of contacts between susceptible and infected individuals by assuming an exponential decay of the ϵ_i terms, viz., (3) Here, τ_i represents the characteristic time-scale describing the increased effectiveness of non-pharmaceutical interventions. These could include restrictions on crowding, improvements in screening procedures as well as increased testing. The net effect of these could be chosen to be different for each of the infectious categories. Although the force-of-infection ultimately involves a product of both the β(t) and the ϵ(t) terms, splitting them out thus provides somewhat more control over the specifics of the non-pharmaceutical interventions as well as the ability to account for new variants that might affect the value of β but not the ϵ-s.

1.4 Estimating R₀

Given the central equations defining our model, we compute the dominant eigenvalue of the next generation matrix [80] obtained from Eq 1. This is the basic reproductive ratio, R₀. This result, derived using a next-generation method outlined in SI, yields (4) In the equation, β is the infectivity parameter. The reciprocal of the λ’s capture the timescales associated with different stages of the disease progression. Together, they determine the serial interval associated with the disease. Note that in the relative intensity of the contacts, the (ϵ)’s are simply constants; likewise, α and μ representing proportions are dimensionless constants. See Table 2 for the description and values of each of these parameter values.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. The transition rates between compartments and the efficiency parameters.

These parameters are fixed during the analysis and their values are taken from [50, 84].

https://doi.org/10.1371/journal.pcbi.1010632.t002

1.5 Estimating R(t)

The time evolution of R(t) is obtained from Eq 2. We obtain the next generation matrix (), incorporating age stratification and contact matrices. Lockdown changes the contact matrices; we take this into account in our estimation of R(t). The computation of the next generation matrix, and thence R(t) corresponding to Eq 2, is discussed in Section 1 in S1 Text. Comparing model predictions with the data generates samples of parameters. From these samples we compute the bounds on R(t) for each region of interest.

1.6 Computing effective R(t) from the data

A standard way to describe infection spread is to calculate the effective reproduction rate R(t) at a given time t. To estimate R(t) we use a Bayesian approach developed by Bettencourt and Ribeiro [81], later modified by K. Systrom [82]. In this approach, given k new cases, the probability distribution of R(t) on a certain day t is: where P(k|R_t) is the likelihood of seeing k new cases given R(t) (this is assumed to follow a Poisson distribution), P(R_t) is the prior, and P(k) is the probability of seeing k cases. The method of deciding the appropriate priors is described in Ref. [81, 82].

1.7 Fixed parameters

Several parameters determine the progress of the pandemic in the INDSCI-SIM model. These include parameters which remain fixed, such as the rate of transitions between the compartments, described above, and a choice for the initial IFR. Other parameters are allowed to vary so as to model the effects of non-pharmaceutical interventions or optimized to fit available data, as discussed below. We assume that asymptomatic patients are less likely to pass on infection to susceptibles by a factor of 2/3, consistent with data from Ref. [67, 83]. The choices of our fixed parameters are provided in Table 2.

1.8 Variable parameters

We maximize a likelihood function to obtain optimal fits to the data, sampling a broad prior distribution in parameter space to assess uncertainties in our description of the disease progression. Our simulations are initiated around 2 weeks prior to the first death being reported. Since the time evolution of any similar compartmental model with constant parameters and no births or deaths would yield only a single wave, we use an adaptive parametrization to address multiple waves of COVID-19 cases. The parameters that enter our description are described below.

We choose a generic functional form for all our parametrizations that involve parameter changes with time. This form interpolates between a high and a low value, and is a function of a characteristic time-scale. We choose a hyperbolic tangent function for concreteness, and because it is easy to specify, but any other suitable function could be used in its place.

1. Initial exposed (E_initial): The exposed people on the first day of simulation, distributed between age groups according to the population fraction.
2. IFR-associated parameters (Δ_IFR, δⁱ(t_initial), δⁱ(t_final) and D_IFR): A number of studies [85–87], suggest that the IFR has decreased over time. This is attributed to an improved clinical handling, new pharmacological treatments such as corticosteroids, non-pharmacological treatments such as proning and simply earlier interventions and finally the potential prophylactic consequences of lower viral load exposure from masking [85]. For concreteness, we choose a specific functional form to describe this decrease, accounting for it by assuming a smooth transition governed by a hyperbolic tangent function. We model the variation of as: (5) where D_IFR is the transition time-point (here, tied to the Hospitalized Fatality Ratio) and Δ_IFR is the characteristic time width of the transition. This introduces four parameters, δⁱ(t_initial), δⁱ(t_final), the timescale Δ_IFR and the transition point itself, D_IFR. Although this formulation is general, allowing for age-bracket dependent variation of the IFRs with time, we choose to allow all age-brackets to have the same functional behaviour, rendering the index i redundant. We take δⁱ(t_final) to be a sixth of δⁱ(t_initial): given an initial IFR of 0.3% this potentially allows the IFR to decay across the range 0.3% to 0.05%, over a timescale defined by Δ_IFR, with the cross-over point between the upper and lower limit defined by D_IFR.
3. Bias(b) and bias-variation (Δ_b): A large fraction of the population remains asymptomatic to COVID-19 infections. In addition, some symptomatic patients may also opt not to be tested. Since testing in India has been limited throughout the first wave, a substantial percentage of actual infections can be expected to have remained undetected, with detected infections always represent an under-counting of the true numbers of infected [88]. We use a bias parameter b to accommodate this scaling relation (actual daily infection = b(t) × reported infection), introducing another parameter to estimate. We parametrise the bias factor as: (6) where Δ_b represents the bias variation timescale. We use the transition time D_bias to be 3 months from the beginning of March but checked that the results should depend minimally on D_bias if the priors on the b(t_initial) and Δ_b are wide enough. This introduces three parameters b(t_initial), b(t_final) and Δ_b. We fix b(t_final) so that it will reach 1 asymptotically. Our calculation obtains the daily number of infected people from the day. We then fix a multiplier that relates the actual numbers of infected to those detected each day. This number, for those expected to be detected in each day, is then compared with the numbers of those reported infected.

If the reported numbers of daily infected cases, or the daily numbers of deaths, shows multiple peaks, a single parametrization will not capture this behaviour. The sources of this behavior are complex functions of government policy constraining the contacts between people, of the sum of individual actions taken to prevent infection and also of the entry of new potentially more infectious variants. Our analysis must be flexible enough to account for these. Set against this is the requirement that we should not over-determine the model, by allowing for a large number of such changes.

We allow our parametrization to vary in the following way: We assume that the time over which we choose to model the data is divided into N segments, with each segment denoted by i. This defines N − 1 break points or nodes. These node points are same for both new cases and deaths, with the lag between these already taken into in the model via various timescales of progression of the disease. The infectivity parameter β and the timescale τ can change in these segments. For simplicity, τ is assumed to be the same between segments. We thus parametrize the model with the following:

1. β_i: Infectivity within the i’th window.
2. τ_i: Timescale within the i’th window.
3. Node_i: The position of the break points (dates) across which β_i and τ_i are change. We decide the required number of break points by comparing the Bayesian evidence for the models with different break points.

Comparisons of the Bayesian evidence indicates that different τ_i’s defined in each segment are not favored. Thus, we work with N β’s, a single τ and N − 1 Node parameters. For N windows, then, we will have N + 1 + (N − 1) parameters specific to adaptive parametrization, 4 parameters (E_initial, D_IFR, Δ_IFR, b, Δ_b) common to both parametrizations and 2 noise error parameters (for modelling the fluctuations in the infected and death data discussed in subsubsection 1.9.1).

1.9 Bayesian modeling

Our model is a mechanistic model that should capture the significant aspects of the dynamics of COVID-19 disease spread. It has a number of parameters which must be optimized. We achieve this optimization through a simultaneous fit to the data on deaths and on detected cases. The latter is related to the “true” number of cases through the time-dependent bias factor.

We address this optimization through Bayesian methods. In such methods, probability distributions over parameters, and not point estimates, are obtained [89, 90]. Bayesian models require a prior distribution over parameters to be specified as well as a likelihood of observing the data under a specific assignment of parameters. Given parameters and initial conditions, a compartmental model defines a unique solution for each of the compartments.

Such methods begin with prior estimates for the parameters entering the model, usually chosen so that they are at least partially constrained by prior knowledge [89]. A likelihood function, chosen appropriately, estimates the probability with which the observed data can be accounted for by a specific parameter set. This leads to a posterior distribution over the parameter set. Credible (confidence) intervals can be derived from such calculations.

We now address the role of initial conditions and sensitivity to parameter values in the context of our Bayesian modelling for parameter estimations. We initiate our model simulations when the infected population is small enough that the susceptible population is not substantially depleted. With this starting point, the dependence on this initial condition is negligible. All other important parameters are estimated using Bayesian methods which incorporate uncertainties in their values. Since we sample over an ensemble of initial values, the dependence on initial conditions is also expected to be minimal. Furthermore, Bayesian modelling naturally allows for parameter variation (e.g. for parameters listed in Table 3), finding an optimum range of posterior distributions constrained by the observed data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Priors on model parameters used in our analysis to account for time varying parameters.

https://doi.org/10.1371/journal.pcbi.1010632.t003

1.10 Data

We use cumulative infected and death data downloaded from the COVID19India website https://www.covid19india.org/. We manually curated the time series, but only wherever required to deal with major outliers arising from reported corrections to the record. (For example Mumbai reported 917 deaths on June 16th 2020, as opposed to the average count of 78.3 for previous week.) These outliers were smoothed by replacing their values by the local average and redistributing the excess count over 20–60 previous days depending upon the value. Typically, for each time series there are only a handful of such outliers (< 5), if at all, and are removed manually. For our simulations, we need the population data for each zone of interest (district, state or country), and mention each of these sources in the appropriate sections. For each such zone, we also need the population fraction in each age group which we obtain from the 2011 census data [95].

1.11 Estimating mortality under-counting

Reported infections are expected to be substantially lower than the actual numbers of infected. While death numbers should nominally be better indicators of the state of the pandemic, most deaths in India happen at home and a medically certified cause of death may not be available, given low levels of MCCD coverage. Moreover, death registration is not uniform across the country. Estimates of overall under-counting of covid deaths vary, but reasonable estimates suggest a factor of 1.5—5 across India during the first wave [49].

Under-counting for cases and deaths cannot be estimated independently in our model, since one of these can be subsumed into a definition of the IFR. Benchmarking our results to those from serological surveys (serosurveys) provides a way of estimating this undercounting, since such surveys provide an estimate of the fraction of the population which has sustained a prior infection by the time of the test. Requiring that the quantum of deaths be consistent with the postulated IFR then allows us to estimate the undercounting of deaths at the time of the survey.

We analyse data for Karnataka—a typical Indian state with both rural and urban districts. (Approximately 62% of Karnataka’s population lives in rural areas, comparable to a figure of 65–70% for India as a whole.) We choose Karnataka for our analysis so that we can compare our numbers to district-wise serosurvey results available from a detailed study conducted in early September of 2020 [96].

For this part of the analysis, we assume that IFR does not change over the course of the first wave; we refer to this as the fixed IFR analysis. We choose several fixed values of (age-averaged) IFR in the range 0.15% and 0.3%. We vary the extent of death undercounting by multiplying the observed number of deaths by a factor >1. We then chose a baseline district where it is potentially safe to assume that death under-counting might be minimal and that, hence, . For this district, we expect the number of infections as estimated from serosurveys should coincide with those obtained from our calculations as well as be consistent with the numbers of deaths given our assumed IFR.

Once we standardize our choice of (age-averaged) IFR with this method for a baseline district, we assume that the same IFR holds true across the state. We can then use the serosurvey results to estimate the true number of deaths for each district [96]. If D_reported is the number of reported deaths, I_sero is the expected number of infected individuals based on serosurvey, we can in principle estimate the death under-counting factor using the relationship between these quantities: . This relation is an average estimation. However, as the pandemic is ongoing during the serosurvey, the infection and deaths differ by a lag (≈ 1–2 weeks) and thus, we cannot use this relation directly.

Instead, using time series of reported infections and deaths, we first estimate actual infections by the time of serosurvey (X(Z) for a district Z). If there is no death under-counting, the expected proportion of actual infections will be equal to the actual proportion of infected from serosurveys (Y(Z)). We estimate any offset from this expectation to estimate a death under-counting factor for each district Z: . All these estimates were performed by analysing time-series data that begins with an initial simulation date (two weeks prior to the available infected and death reports of each district) and goes on to early September 2020, when the Karnataka serosurvey was conducted. We simultaneously fit the daily time-series of infection and deaths, using the fixed values of the (age-averaged) IFR mentioned above. Uncertainties in the serosurvey results translate to uncertainties in our determination of . Apart from the estimation of death undercounting, the rest of our analysis is based on a time-varying IFR described in the previous sections.

We also employ a ‘death multiplier’ for our analysis of Karnataka state, which is different from . The death multiplier is provided as an input that represents our estimate the amount by which deaths are under-counted. We multiply the death numbers uniformly by the death multiplier for the entire period of analysis. For consistency, these numbers must be aligned to the IFR and the underlying numbers of infections. is the parameter that we obtain as an output that indicates possible death undercounting. For example, in a particular case, if without the use of death multiplier, we find , then it indicates that the death numbers are under-counted by a factor of 2. For the same case, if the death multiplier was input was 2, then we will estimate .

2.1 Karnataka districts: Estimating death under-counting

Our analysis of the Karnataka data follows a relatively complex set of protocols, which we use to obtain reasonable estimates of death undercounting and total infected in the first wave. This is a challenging exercise because (a) the cases are under-counted, (b) deaths are under-counted and (c) there is no true consensus on the IFR values that should be appropriate. Our method deals with these issues by assuming that the deaths are calibrated correctly in one urban district with the best health facilities in the state.

We begin with an estimation of death under-counting in Karnataka (subsubsection 2.1.1), using a simultaneous fitting of reported infections and deaths to our INDSCI-SIM model. We then compare this with serosurvey findings, following the fixed IFR analysis protocol described in subsection 1.11. Since the serosurvey we use to calibrate was conducted in September 2020, this first analysis is carried out for daily reported data of new cases and deaths from the start of the pandemic in India to early September 2020. We find that an IFR of 0.2% seems a reasonable estimate upto early September 2022. We then apply our model allowing for time-varying IFR (and other parameters, as described in subsection 1.1 to 1.10), to data of each district of Karnataka (subsubsection 2.1.2); we find that different districts had different levels of infections, ranging from 20% to 70% in some cases, with a state-wide aggregate under-counting factor of 2.2. Finally, we present results for the aggregate state-level data for Karnataka (subsection 2.2).

2.1.1 Estimates of death under-counting upto early September 2020, with fixed IFR.

We use the serosurvey results published in Ref. [97] for all districts of Karnataka. We first scale all age-dependent IFR’s with a constant factor to ensure age-averaged IFRs of 0.3%, 0.2% and 0.15%—these have been argued to be in the right range for LMICs [78]. We assume further that either (a) there is no death under-counting or, (b) that deaths have been under-counted by a factor of two. As can be seen in Table 4, assuming that deaths have been under-counted by a factor of two doubles our estimate of the actual infected in the same IFR category, as expected. On the other hand a decrease in the IFR results in an increase in our estimate of the actual numbers of infected. We stress that, for the purposed of this specific analysis, the IFR is assumed to be constant over time (no change in δⁱ) and that we only analyze the data till September 15, 2020 for all districts in Karnataka for consistency with the Karnataka serosurvey. Population data for Karnataka districts is sourced from the projections for the period 2020–2021, as provided by the report issued by Directorate of Economics and Statistics, Bangalore [98]. We choose Bengaluru Urban as our baseline model district, assuming that the recorded numbers of deaths are accurate. This choice is motivated by the observation that Bengaluru Urban is the largest district in Karnataka by population. It contains the state capital and is also thus likely under greater scrutiny than other, more remote, districts. Comparing our model results with those from the serosurveys for Bengaluru Urban we find that an IFR of 0.2% (highlighted column in tha table) yields numbers of those actually infected that are close to the measured seroprevalence (the undercounting factor ). To the extent that the reported death numbers are accurate for Bengaluru Urban, this is then our best estimate for the IFR—but only up to the date of the serosurvey. We can now use this assumption to estimate levels of undercounting in other districts, . (We, however, expect that this IFR value will go down further as time goes on, reflecting improvement in patient treatments and a consequent mortality reduction as well as COVID deaths in older age-bands leading to a reduction in the population which carries the most mortality risk.)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Estimation of death under-counting by comparing seroprevalence.

The first column lists each district of Karnataka (as on Sept 2020). The second column provides the seropositivity estimate as obtained in a serosurvey across the first two weeks of September 2020. The third column provides our estimate of the potential death undercounting in that district, , represented by the ratio of seropositivity to the actual number of reported infected. Since death overcounting is not possible, the minimum value of can only be 1. The uncertainty on this number is obtained by propagating the 95% uncertainties from the serosurvey reports to our prediction of actual infection. In cases the lower bound on goes below 1, we present the mean value and the range from 1 as an upper bound in parentheses. In a few cases we just present 1 as an upper bound, representing those cases where even the upper limit is less than or equal to 1.

https://doi.org/10.1371/journal.pcbi.1010632.t004

Using these columns in Table 4, for IFR = 0.2%, we plot the factors in Fig 3 as a choropleth in the left panel. Darker colors for districts represent a larger death undercounting. In the right panel we plot for IFR = 0.3% analysis. As expected, the estimated death undercounting increases here. This would be consistent with an death undercounting in the Bengaluru Urban district by 60% i.e. ().

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Ratio of seroprevalence reported in districts of Karnataka across the first two weeks of September, 2020 and our prediction of actual infected on the same date.

This ratio indicates the death under-counting in different districts assuming uniform IFR in different districts. We plot the ratios here for IFR = 0.2 [a: left plot] and 0.3 [b: right plot] respectively. Note that these numbers are plotted from Table 4 for corresponding IFRs. For IFR = 0.3 the INDSCI-SIM prediction of actual infected is certainly lesser than the IFR = 0.2 case. The plot on the right shows larger ratios compared to the left. The maps are generated using the data made available at https://github.com/covid19india/covid19india-react/tree/master/public/maps. As per their declaration the repository is available under an MIT License: https://github.com/covid19india/covid19india-react/blob/master/LICENSE.md.

https://doi.org/10.1371/journal.pcbi.1010632.g003

2.1.2 Estimates of cases by the end of first wave (Feb 2021), with variable IFR.

Having estimated , we redo our analysis for each of Karnataka’s districts, multiplying reported deaths by the multiplying factors tabulated in Table 4. We then run the program forward till February 15, 2021, but now also allow the IFR to vary according to the parametrization described earlier. Using this analysis we tabulate the actual infected percentages we predict by Dec 15, 2020 and Feb 15, 2021 in Table 5, together with 95% confidence intervals. We find that between 20% and 70% of the population in these districts has been infected by February 15, 2021. We find no substantial change in the actual infection between December 2020 and February 2021.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Estimation of actual infection in Karnataka districts on December 15, 2020 and February 15, 2021.

The uncertainties are quoted at 95% C.L. This table is obtained from an analysis that assumes a time varying IFR as explained in Methods.

https://doi.org/10.1371/journal.pcbi.1010632.t005

The relationship between actual numbers of prior infections and results from serosurveys is complicated by the observation that antibodies measured in such surveys have been observed to decay in concentration over a period of several months. However, given the fact that the peak in cases in Karnataka occurred over September-October 2020, and that the serosurvey was performed during this period, it is reasonable to assume that this will not impact our interpretation of the serosurvey observations [97].

2.2 Projections and data fitting for Karnataka state: The effect of death undercounting

We now present our analysis of Karnataka state viewed as an aggregate of all its districts. Karnataka has both urban and rural districts. We will assume that a suitable estimate of undercounting at the level of the entire state can be obtained through a population average of the undercounting at the level of its individual districts. Accounting for the fraction of the total population associated to each district, multiplying them by the district undercounting factor U(D) and summing these values gives us a population averaged death undercounting for Karnataka to be 2.2.

In this section we analyze the Karnataka infection curve in two ways. In the first, we assume no death undercounting. In the second we multiply the reported deaths by 2.2 to account for state-wide undercounting. Here the reported infection and deaths data are addressed using our model. We have used the data till the end of the first wave, i.e. February 15, 2021. From the sampled parameter space we predict the time evolution of bias (between actual and reported infection), the IFR and the actual infected population.

We present the results via a set of 6 plots in three panels. In the top panel we present the fit to the daily reported cases (left—(a)) and to the daily reported deaths (right—(b)). The actual predicted infections and the bias are plotted in the middle left (c) and right (d) plots respectively. Note that the bias factor, once multiplied with the reported infections, provides the estimate of actual infection. The change in the IFR (left—(e)) and the effective reproduction number R(t) (right—(f)) estimates are plotted in bottom panel. For R(t) we also plot an independent estimate [82] as a reference.

2.2.1 Karnataka: No death undercounting.

In Fig 4 we provide model results for Karnataka state data assuming deaths have been counted accurately. The posterior distributions and parameter constraints are provided in Fig F and Table A in S1 Text. Both infected and death numbers are well described by the model. The state-wide data has a single peak in both infected and deaths, indicating that the many variable peaks present in the individual districts are averaged out to yield the relatively simple structure shown in the figure. The numbers for those totally infected ranges between 15–45% at the 2σ level by mid February. We see that the strong urban to rural gradient in seropositivity, once averaged across the state, yields infected fractions that are below those in the urban regions. The bias multiplier follows results obtained from the city of Bengaluru.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Timeseries analysis for Karnataka without death undercounting.

We plot the fit to the daily infected cases [a: top left] and daily reported deaths [b: top right] assuming no death undercounting. Middle panel contains the cumulative actual infected cases [c: left] and the bias multiplicative factor [d: right] obtained as a ratio between actual and reported infections. The left plot [e] at the bottom panel contains the evolution of age averaged IFR. The bottom right plot [f] contains our estimation of R(t) and an independent [82] measurement. Note that the bands correspond to 2σ and 3σ confidence levels. We do not assume any death undercounting in this analysis.

https://doi.org/10.1371/journal.pcbi.1010632.g004

The IFR shows a decreasing trend after lockdown with a final range of IFR between 0.05–0.12%, consistent with estimates from the serological surveys [97]. These surveys find an overall adjusted total prevalence of 27.7% (95%CI 26.1–29.3); we obtain about 25% at the mid-line of the 95% confidence interval in mid-September. The case-to-infection ratio was estimated to be 1:40 while the infection fatality rate was 0.05%; our results are consistent with both observations. The R(t) band shows an increase after the lockdown till August 2020. After that R(t) decreases till December 2020 and stays below 1 till February 2021. In this case, bounds on cumulative infections and deaths are plotted with data in Figs C-E in S1 Text.

2.2.2 Karnataka: Incorporating death undercounting.

Fig 5 shows our model predictions for the Karnataka state data where the death timeseries is uniformly scaled with multiplier of 2.2. The related shifts in the posterior distributions and parameter constraints can be seen in Fig F and Table A in S1 Text). The daily infection and the scaled daily death report fits are similar to Fig 4. However, compared to the analysis that assumed no death multiplier, predictions for the fraction of those infected ranges between 30–50% Table 5. A larger death undercounting thus leads to a higher estimates of actual infection within the population. Correspondingly, the bias multiplier is also large.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Timeseries analysis for Karnataka with death undercounting.

We plot the fit to the daily infected cases [a: top left] and daily reported deaths [b: top right] taking into account estimated death undercounting. Middle panel contains the cumulative actual infected cases [c: left] and the bias multiplicative factor [d: right] obtained as a ratio between actual and reported infections. The left plot [e] at the bottom panel contains the evolution of age averaged IFR. The bottom right plot [f] contains our estimation of R(t) and an independent [82] measurement. Note that the bands correspond to 2σ and 3σ confidence levels. In this analysis the reported death numbers are multiplied with 2.2, the average death undercounting we obtained from Karnataka districts assuming Bengaluru Urban does not have any death undercounting.

https://doi.org/10.1371/journal.pcbi.1010632.g005

With higher death numbers, the effective IFR decreases at a later stage (around September 2020), converging around February to the values without death undercounting. The effect of death undercounting is absorbed mainly by the bias multiplier and therefore is reflected in the actual infected numbers. We find that the R(t) is nearly identical to the analysis without death undercounting. Assuming death undercounting, our estimates for the cumulative infections and deaths are plotted alongside the data in Section 5, Figs C-E in S1 Text.

2.3 Infection curves for select Indian cities

In this sub-section, we analyse of the progress of COVID-19 across the period Jan 2020—February 15, 2021, for a few select Indian cities: Bengaluru, Chennai, Delhi, Pune and Mumbai. For each of these cities, the data correspond to the urban district agglomeration within which the city is embedded. We assume that deaths have been counted largely accurately for these cities, and thus that U(D) defined above is 1. We will return to the implications of this assumption later.

All these cities show more than a single peak in reported cases. The origins of these peaks are complex. Together, they represent a combination of relaxations of the lockdown and of physical distancing, festivals and political activity at various times, the ebb and flow of migration from nearby districts, as well as day-to-day variations in testing. These effects are captured in our multiple window model with adaptive selection, since infectivity parameters can now vary across different windows and the bias factor between actual and reported cases can itself evolve.

We present a set of 6 plots in three panels, as in results of Karnataka in the previous subsection. We show the results for Mumbai in a subsection below, and relegate the rest to Tables D-G and Figs K-V in S1 Text.

2.3.1 Mumbai.

Results for the city of Mumbai are provided in Fig 6. The population of the district is taken from the projections for 2021 [99]. The posterior distributions and related parameter constraints are provided in Fig H and Table B in S1 Text. The data for Mumbai show two well resolved peaks around June-July 2020 and September-October 2020. Following the second peak, reported infections reduce and remain roughly constant from about December 2020.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Timeseries analysis for Mumbai.

We plot the fit to the daily infected cases [a: top left] and daily reported deaths [b: top right] assuming no death undercounting. Middle panel contains the cumulative actual infected cases [c: left] and the bias multiplicative factor [d: right] obtained as a ratio between actual and reported infections. The left plot [e] at the bottom panel contains the evolution of age averaged IFR. The bottom right plot [f] contains our estimation of R(t) and an independent [82] measurement. Note that the bands correspond to 2σ and 3σ confidence levels.

https://doi.org/10.1371/journal.pcbi.1010632.g006

The reported deaths shows a peak and a plateau before reaching a second, smaller peak. While the second peak in the reported infection is higher than the first peak, the plateau in the reported death numbers around October 2020 is substantially lower than the July 2020 peak, indicating better treatment interventions. There are no systematic outliers in the death numbers, once we apply our smoothing to the data. The sharp decrease in death cases in contrast to the increase in reported infections could result from increase in testing (a lowering of the TPR) and/or a decrease in IFR.

While the bias multipler decreases from 100 to 15, the IFR does not decrease as rapidly. At 2σ the data is consistent at about an IFR of 0.15%, which is the largest when compared to the previously discussed cities. In Mumbai, a relatively smaller fraction of population (50%) appears to be infected within this calculation. However, fluctuations in the infected data w.r.t. the mean value translate to the uncertainties in the time-series estimation. Thus, the bands are wider compared to Bengaluru, Chennai and Delhi. At the 2σ-level we estimate that 54–84% of total Mumbai population may have been infected. R(t) shows a steady decrease within the lockdown and after that till mid-August 2020 when it starts to rise again temporarily indicating an upcoming second peak. Following the peak R(t) first reduces below 1 for a few months, then rises to fluctuate around 1 in early 2021 till February.

Mumbai serosurvey data showed that 54·1% of samples in slums and 16·1% of those in non-slums tested positive [100], in a study performed between June 29 and July 19, 2020. These are broadly consistent with our estimates. Our estimates for cumulative infections and deaths are plotted with data in Fig G in S1 Text.

2.4 India

To describe the data for India, we assume that the death undercounting factor that represents India as a whole is the same as that obtained for the state of Karnataka, a factor of 2.2. In doing this we are motivated by the observations that the urban-rural gradient can be substantial and that the division between urban and rural populations at the Karnataka state level and at the national level are not too different, at an approximate ratio of about 35:65.

The national data for COVID-19 infections shows a single peak structure around mid-September 2020, at a bit less than 98,000 cases. The daily numbers then began to reduce after that reaching a value of about 10,000 by January of 2021. Our model fits infection and death reports in three adaptive windows between March 1, 2020 and February 15, 2021. For our calculations, we use India’s population data from the projections reported in the 2019 report published by National Commission on Population [101].

We find the bias multipliers change from 120 to 20 in this time period, thus indicating that by the end of the first wave, infections were being underestimated by the reported cases through a factor of 20. This decrease is in accordance with the reported overall increase in testing during this period. At the 2σ level, by the middle of February when the second wave began, we find that between 20–50% (95% CI) of Indian population was infected, with a point estimate of about 35%. This is a more conservative estimate than made by others. However, we believe it is consistent with the third ICMR serosurvey across December and January as well as with the slow pace of case accumulation over January and February.

The IFR band is consistent with a global India IFR of about 0.1% (95%CI 0.05—0.15). The R(t) band has several interesting features corroborated with the independent measurement plotted in dots. Before the lockdown, the R(t) stays in the high range 2.5–3. During lockdown it reduced to 1.5 and thereafter to the vicinity of 1.2. It stayed around 1.2 till August, going below 1 around October, signaling a steady decrease in infected numbers. A slight upward blip in November 2020 is consistent with the festival season in the north and east of India. The R(t) then stayed below 1 through till February 2021. Assuming death undercounting, estimations for cumulative infections and deaths for India are plotted with data in Fig I in S1 Text. We present the posterior distribution of model parameters for the national data analysis in Fig J in S1 Text.

Our estimates indicate that at the beginning of March the mean exposed population in India was around 6000. Around 14 March the number of officially reported infections in India crossed 100, but we expect that a lot of infections must have escaped detection.

Three adaptive windows are defined between simulation onset (March beginning) to 154’th day from the onset (August 2, 2020); between August 2, 2020 to October 21, 2020 (234th day from the simulation onset) and October 22, 2020 to February 15, 2021 (simulation end). These 3 regions have three different associated infectivities β₁ = 0.0755, β₂ = 0.117 and β₃ = 0.220. This is the base infectivity with respect to the simulation onset. This is suggestive of two phases in which increased relaxations of COVID-19-associated restrictions may have led to a rise of cases, one around the beginning of August and one around the end of October, although one should be careful not to overinterpret this data,

Note that the effective β is modified by the parameter τ which has a mean value of 142 days. Therefore while the β is increasing with different windows, the effective β is also increasing. The bias parameter is unbounded from above even with a conservative prior. However the time evolution of bias determined by b and Δ_b for national data indicates a steady decrease, reaching a bias of around 20 as plotted in Fig 7. The IFR shows a transition around mid July (108 day since the onset) with a width Δ_IFR of 2 months. σ₁ and σ₂ represents the noise in the infected and deaths data respectively. They are well constrained.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Timeseries analysis for Indian national data.

We plot the fit to the daily infected cases [a: top left] and daily reported deaths [b: top right] assuming taking into account estimated death undercounting. Middle panel contains the cumulative actual infected cases [c: left] and the bias multiplicative factor [d: right] obtained as a ratio between actual and reported infections. The left plot [e] at the bottom panel contains the evolution of age averaged IFR. The bottom right plot [f] contains our estimation of R(t) and an independent [82] measurement. Note that the bands correspond to 2σ and 3σ confidence levels. In this analysis the reported death numbers are multiplied by 2.2 to take into account a measure of presumed undercounting.

https://doi.org/10.1371/journal.pcbi.1010632.g007

The triangle plot of posteriors reveal strong correlations between parameters. Given the time-series of infected persons, a fit to the initial data can be addressed either by having a large initial exposed and lower infectivity or by having a smaller exposed population and a higher infectivity. Such a negative correlation is seen in Fig J in S1 Text, between E_initial and β₁. Since smaller τ₁ results in a faster decrease in infection by lowering the effective β, to fit the infection data, if β_i increases, τ must decrease to maintain a similar level of effective infectivity. These negative correlations are reflected in the contour plots as well. The slope of the infection curve is given by β; the β’s are thus positively correlated.