Persistence of SARS-CoV-2 immunity, Omicron’s footprints, and projections of epidemic resurgences in South African population cohorts

3.1 Seroconversion and vaccination

Combining the longitudinal rRT-PCR assays and serological test results from the two study sites, we found that the seroconversion rate was 97.6% (583/600) among primary infection episodes with at least one blood specimen collected 30 days after their first rRT-PCR positive test. Among 639 individuals with negative serological specimen(s) after their first blood specimen was taken but who seroconverted later, 86% (549/639) were confirmed by rRT-PCR during the study period, suggesting that twice-weekly rRT-PCR testing captured most of the infections during the study period. Among 706 seroconverted individuals with at least one follow-up blood specimens after seroconversion, only 40 (5.7%) later sero-reverted (sero-reversion is defined as sero-converted individual who became sero-negative the subsequent serologic test). The baseline characteristics of all individuals are reported in Table 1. 10% individuals received at least one dose of SARS-CoV-2 vaccine and 5% were fully vaccinated within the study period of PHIRST-C.

3.2 Defining and typing variant-specific infection episodes

For each individual in the study, we defined the duration of a SARS-CoV-2 infection episode as the time interval between the first and the last of a set of positive rRT-PCR test, where consecutive positive tests were separated by less than 30 days. If at least one positive rRT-PCR test within the infection episodes was identified as any of the variants of concern (VOC) (defined by the WHO) Alpha (B.1.1.7), Beta (B.1.351), Delta (B.1.617.2), the lineage of the infection episode was assigned to the identified VOC. If all lineage-typed rRT-PCR tests within an infection episode were identified as the D614G, then the lineage of the infection episode was assigned to D614G. If none of the positive rRT-PCR specimens within an infection episode had a lineage defined, the lineage of the infection episode was labeled as inconclusive. In total, we observed 669 infection episodes, including 634/669 (95%) with defined lineages. For the 35/669 (5%) infection episodes with inconclusive lineages, we assigned the infection to the dominant SARS-CoV-2 lineage at the study site identified within the month of the infection episode. After lineage extrapolation, 108/669 (16%) of the infection episodes were estimated to be caused by D614G lineage; 245/669 (37%) infection episodes were estimated to be caused by VOC Beta lineage; 299/669 (45%) infection episodes were estimated to be caused by VOC Delta lineage; 17/669 (3%) were estimated to be caused by other lineages.

3.3 Characterizing the relative viral RNA shedding kinetics of SARS-CoV-2 infection episodes

When interpreting Ct values from rRT-PCR test by the Seegene Allplex™ 2019-nCoV kit, it is important to stress that the Ct value of a single rRT-PCR test is not a direct measurement of the quantity of viral genetic material present in an individual specimen (in absolute terms) (57). Many factors could influence the Ct value of a rRT-PCR test, including but not limited to the specimen quality, extraction method, chemistry of reagents, gene targets. Further, Ct values cannot be directly compared between assays of different types (58). However, comparing serial Ct values and/or Ct values of different groups of population collectively, based on rRT-PCR tests from the same assay and laboratory setting, does reflect the relative variation in terms of viral genetic material concentration over time and between population subgroups (57, 58).

We measure the viral RNA shedding intensity of a given specimen as the cycle threshold value of Seegene Allplex™ 2019-nCoV kit’s N gene target Ct value of the specimen, i.e., . For the Ct value, we are only considering the N gene target to avoid between-target variation. In this study, instead of directly interpreting the Ct value of a single rRT-PCR test, we fit a mathematical model to capture the temporal kinetics of Ct for longitudinally collected specimens and extract statistical summaries of RNA shedding intensity for each infection episode. In particular, following the method presented in (21), for each infection episode, we modelled the Ct kinetics during the “RNA proliferation stage”, characterized by linear decrease in Ct until a trough in Ct is reached; and during the “RNA clearance stage”, characterized by linear increase (in Ct) from the trough until the last rRT-PCR positive result. The duration of the “RNA proliferation stage” τ_p was defined as the time from Ct first exceeding the detection threshold (i.e.,) until the trough of Ct; the duration of the “RNA clearance stage” τ_c was defined as the time from trough Ct to Ct reaching above detection threshold (i.e.,). The duration of rRT-PCR positivity τ_s for each infection episode was defined as the total duration of the RNA proliferation and RNA clearance stages: τ_s = τ_p + τ_c. The duration of RNA proliferation τ_p, RNA clearance τ_c, rRT-PCR positivity τ_s, as well as the trough shedding intensity were estimated for each infection episodes using Markov Chain Monte Carlo (MCMC) method.

The fit of the RNA shedding trajectories for each variant are shown in Figure 2A-G. In addition, we performed multivariable regression analysis to evaluate the dependencies of duration of rRT-PCR positivity and trough Ct on participants’ characteristics including age, sex, BMI, HIV infection status, symptom status, variant type, and evidence of prior infection. Since the nasal swab sampling period ended on August 28, 2021, around the peak of the Delta wave in both sites, we limit our analysis to infection episodes with first positive PCR sample 30 days prior to the end of sampling to avoid censoring bias. The result of the regression is presented in Figure 3A-B.

3.4 Assigning the lineage of prior infections among seropositive individuals

For individuals who seroconverted without a rRT-PCR confirmed infection episode (prior to the start of the PHIRST-C cohort), we assigned the lineage type of the individual’s unobserved infection according to lineages’ prevalence (based on the infection episodes that had lineage information) at the study site in the month of the earliest seropositive specimen.

3.5 Risk factors associated with SARS-CoV-2 (re)infections

For each individual in the urban and rural sites, we first reconstructed the status of prior and ongoing infections at a daily time resolution, based on the lineage-typed infection episodes and serologic results described in Section 3.3 and 3.4 above. Specifically, we denoted a person-day observation of individual i’s infection status at day t as. We set if individual i acquired infection on day t (i.e., the date of viral acquisition, marking the start of an infection episode), and set if individual i didn’t acquire infection on day t. Individuals without rRT-PCR positive results throughout the study period were assumed to escape infection throughout the observation period (, for all t). For individuals with observed infection episodes during the study period, the individual Ct kinetics were modelled as in Section 3.3. The model allowed us to estimate the time at which Ct crossed the rRT-PCR detection limit (Ct = 40). We assumed that the time t_inf of acquiring infection occurred δT_cryptic = 4 days prior to the Ct crossing the detection threshold (59). The time period from infection onset to 30 days after viral RNA clearance (i.e., the estimated time of Ct crossing below the rRT-PCR detection limit at cycle threshold of 40) was considered a period of active infection, when an individual cannot get infected again. Thus, this refractory period was excluded from the survival analysis of (re)infection hazard. We also considered a sensitivity analysis with a 15-day exclusion after viral RNA clearance. We also censored participants after they received their first dose of vaccine. In Figure S1, we visualize the vaccine uptake overtime for both J&J/Janssen Ad26.COV2.S and the Pfizer/BioNTech BNT162b2.

After having reconstructed the daily infection status of each individual, we modeled the risk of (re)infection for each individual using piecewise exponential models of survival analysis with time varying covariates (60, 61). The piecewise exponential model is estimated by performing Poisson regression with the daily infection status of each individual as the binary outcome and daily covariates (60, 61). The covariates considered in the regression analysis included the individual’s age (allowing for variant-specific age effects), sex, body mass index (BMI), HIV infection status, household size, household crowding, variant type, study site, time since prior infection, SARS-CoV-2 exposure intensity from household members with on-going SARS-CoV-2 infection, and SARS-CoV-2 prevalence in the community. We estimated that the SARS-CoV-2 exposure intensity from household members with on-going SARS-CoV-2 infections at time t as the sum of the estimated dCts of all infected household members (62). The dCt of household member j at time t was measured as . Thus, the household exposure intensity exerted on individual i at time t can be expressed as , where j is the sum over all infected household members at time t. We further discriminated the household exposure intensity during the RNA proliferation and RNA clearance phases of infections in household members. We estimated the weekly SARS-CoV-2 community prevalence as the prevalence of each variant in the entire study site (no of variant-specific infections that week/population tested that week) and use the variant-specific prevalence as a proxy for the SARS-CoV-2 community exposure intensity. We included household- and individual-level hierarchical random effects, as well as “day-of-week” and “day-of-year” random effects in the regression model. The full model selection procedure is reported in Table S1 (the model selected was model 16, corresponding to regression results presented in Figure 3C).

We considered whether saturation effects could possibly affect our estimates of infection risk in the proliferation and clearance phases. If the risk of transmission in the household setting was very high during the proliferation stage, so that most household members were infected at the end the proliferation stage of the index case, there may not be enough observations to power the estimate of transmission hazard during the clearance stage. Empirically, however, we found that transmission saturation was uncommon for the participating households. We define a household transmission cluster as a group of infections separated no more than 14 days by their time of infections (equivalent to hierarchical clustering based on infected individuals’ time of infection with single-linkage). In Figure S3, we plot the distribution of the size of household clusters as a function of the number of individuals in the household (household size). We found than only 10 household transmission clusters were saturated (all members were infected), and the median size of the household cluster was 38% (proportion of household members infected). Thus, saturation was not an issue, and there were enough susceptible individuals to power the estimation of the hazard of transmission form infected individual during the clearance phase of infection.

Since the nasal swab sampling period ended on August 28, 2021, around the peak of the Delta wave in both sites, we limit our analysis to infection episodes with first positive PCR sample 30 days prior to the end of sampling to avoid censoring bias. In total, 21% (303407 person-days/1472400 person days) of the total person-days of observation were excluded from the regression during the entire study period due to missing nasal swab visits, missing serologic status, or individuals experiencing an active infection episode. A household with an individual with chronical SARS-CoV-2 infection was also excluded. The hazard ration of the piecewise exponential model is presented in Figure 3C (corresponding to Model 16 in Table S1). We conducted two additional sensitive analyses with slight variation of this model: in the first sensitivity analysis (Figure S4A) we censored individuals until 15 days post viral RNA clearance (instead of 30 days in the main analysis); in the second one (Figure S4B) we used logistic regression (logic link) rather than Poisson regression (log-link) (63).

4.1 Estimating the reporting rates of D614G, Beta, and Delta during the first three epidemic waves in the Dr Kenneth Kaunda District

Here we compare the carefully monitored weekly PHIRST-C infection data with weekly surveillance data for the broader district to estimate aspects of SARS-CoV-2 dynamics in the district, including the rate of under-reporting for each of the pre-Omicron variants, pre-Omicron infection histories, and the impact of the Omicron wave. The dark blue bars in Figure 4F top panel shows the weekly incidence I_case(t) of SARS-CoV-2 cases in the Dr Kenneth Kaunda District, North West Province, South Africa during the first three epidemic waves captured by SARS-CoV-2 surveillance. We calculate the weekly cumulative case rate C_case(t) at time t as , where week 12, 20202 is the week when the district starts reporting SARS-CoV-2 cases. We denote the cumulative infection attack rate at time t (proportion of the population infected by each of these variants at time t) for D614G, Beta, Delta (and others) variant as , and respectively. We estimate , and within the time period of the PHIRST-C cohort based on PHIRST-C’s serology and variant-typed SARS-CoV-2 infection episodes (Detailed in Materials and Methods Section 3). We denote the reporting rate of D614G, Beta, and Delta (equal to Others) in the urban site as β_D614G, β_Beta, and β_{Delta & Others}. We assume that the reporting rate for each variant (within each epidemic wave) is constant and that Other variants (low frequency, mostly in between second and the third waves. See Figure 1F) have the same reporting rate as Delta. The cumulative case rate C_case(t) and cumulative infection rates , and at time t satisfy . Given , and for t within the PHIRST-C cohort, we can estimate β_D614G, β_Beta and β_{Delta & Others} using linear regression without intercept. The regression analysis was performed using R package lme4 version 1.1-27.1 (64).The estimates are reported in Table S2.

4.2 Estimating the weekly incidence rate of SARS-CoV-2 infection based on case incidence and reporting rate in Dr Kenneth Kaunda

Based on the variant-specific reporting rate β_X (X ∈ {D614G, Beta, Delta & Others}) estimated in Section 4.1 and the variants’ proportions p_x(t) at a given time point t, we can express the overall reporting rate β_Overall (t) at time t as:

Where I_case(t) is the weekly incidence rate of SARS-CoV-2 cases reported to the Dr Kenneth Kaunda District at time t and I_infection (t) is the weekly incidence rate of SARS-CoV-2 infections. We can thus estimate I_infection (t) as

The estimated I_infection (t) prior to the end of PHIRST-C (September 2021) is visualized in Figure 4F top panel (light blue bars prior to September 2021).

4.3 Reconstructing the SARS-CoV-2 antigen exposure history in Dr Kenneth Kaunda by September 2021

In addition to allowing to trace SARS-CoV-2 infection history in detail, the PHIRST-C cohort also recorded participants’ timing of vaccinations. Specifically, at time t, we denote the proportion of the PHIRST-C urban site population with a single prior SARS-CoV-2 infection, which is to D614G, as with and representing the same quantities for Beta and Delta. We denote past vaccination (at least one dose) as with repeat exposures (including repeat infections, vaccination followed by infection, or infection followed by vaccination) as . Then we can express the proportion of population with past SARS-CoV-2 infections at time t as:

If we denote the cumulative infection attack rate in Dr Kenneth Kaunda and assume the PHIRST-C urban site is a representative survey of the population in Dr Kenneth Kaunda, we can express the cumulative prevalence of past exposure of a given exposure type X at time t as:

Where X ∈ { D614G, Beta, Delta, Repeat Exposures, Vaccinated }. In Figure S6 (below) and Figure 4F bottom panel, we visualize the proportion of population with a specific SARS-CoV-2 antigen exposure history.

5.1 Overview of protective immunity

To project the impact of Omicron, we need to model how immunity from infection with pre-Omicron variants and/or vaccination will impact the probabilities of infection, transmission and severe outcomes. The effectiveness of protective immunity (IE) can be measured very broadly as IE = 1 − RR where RR is the relative risk of an outcome of interest (infection, transmission, hospitalization, etc.) and the comparison is against individuals with no prior immunity or with different types of immunity. In particular, we consider three aspects of protection from prior infection/vaccination:

1. Prior infection/vaccination could reduce the host’s susceptibility to reinfection IEⁱ, measured as IEⁱ = 1 − RRⁱ where RRⁱ is the relative risk of reinfection/breakthrough infection when compared to the infection risk in a naïve population, controlling for same level of exposure.

2. Prior infection/vaccination could reduce the risk of onward transmission given reinfection/breakthrough infection IE^{t |i}, measured as IE^{t |i} = 1 − RR^{t |i} where RR^{t |i} is the relative risk of onward transmission for reinfections/breakthrough infection when compared to that of primary infections, conditional on the same contact rate. Here RR ^{t |i} = RR ^td|i × R ^ts|i can be further broken down into the product of reduction in the duration and intensity of shedding (RR ^td|i and RR ^ts|i).

3. Prior infection/vaccination could reduce the risk of disease given reinfection/breakthrough infection. The concept of COVID-19 disease is generic and could encompass a wide spectrum of severity endpoint including symptomatic cases, hospitalizations, and deaths. In this study, we used symptomatic illness as the severity end point. The effectiveness of protective immunity against being symptomatic case, conditional on infection can be measured as IE^{c |i} = 1 − RR ^{c |i}, where RR ^{c |i} is the relative risk of disease for reinfections/breakthrough infections compared to primary infections.

5.2 Modelling boosting and evasion of protective immunity in the Omicron era

Repeat exposures to SARS-CoV-2 antigens can result from repeated infections, booster shots following primary vaccine schedule(s), infection following vaccination or vice versa. Multiple exposures could stimulate a recall response and boost the level of protective immune effectiveness IE through reducing RR further from the baseline protective immunity provided by primary exposure. For example, let’s denote the protective immunity conferred by primary infection with a pre-Omicron (pOm) strain as IE_pOm (pOm) = 1 − RR _pOm (pOm). We can express the level of protective immunity conferred by reinfection with pre-Omicron strains as , with ρ_pOm > 0 measuring the degree of immunity boosted by reinfection with a pre-Omicron strain. For RR_pOm (pOm) < 1), the larger ρ_pOm, the greater the reduction RR_pOm (pOm),when compared to RR_pOm(pOm), and the higher the protection conferred by boosting.

On the other hand, in face of an antigenically distinct variant like Omicron, prior immunity may not be as efficacious due to reduced ability to recognize the antigen through immune memory, leading to elevation in RR and consequently reduction in IE. For example if we denote the protective immunity conferred by primary infection with a pre-Omicron strain against refection by pre-Omicron as , we can express the level of protective immunity conferred by pre-Omicron strain infection against Omicron (Om) as , where σ_Om measures the degree of immune evasion by Omicron. When σ_Om = IE_pOm(Om) = IE_pOm(pOm), indicating no immune escape. When σ_Om = 1, IE_pOm(Om) = 0, indicating 100% immune escape. We can further model the combined effects of boosting and immune escape. For example, if we consider an individual first infected with a pre-Omicron strain, followed by an Omicron infection, we can express the individual’s level of protective immunity against further Omicron infection as .

5.3 Notation conventions and an exhaustive list of different types of protective immunity considered

Following Section 5.1 and 5.2, we use 1 − RR as a measurement of the level of protective immunity. In particular, can be broken down into three independent dimensions: 1) the type of protective immunity - (superscript), 2) the antigen exposure history Y that confers immunity (subscript), 3) the viral strain X against which immunity is directed (brackets).

Where X ∈ {i, t|i, td|i, ts |i, h|i} and:

• i denotes protection in terms of susceptibility against infection.

• t|i denotes protection against transmission given infection.

• td|I denotes protection against transmission through reduced duration of shedding, given infection.

• h|i denotes protection against transmission through reduced intensity of shedding, given infection.

• h|i denotes protection against hospital admission given infection.

Where Y ∈ {pOm, Om, pOm − pOm, − Om } and:

• pOm denotes immune histories including only one antigen exposure by either pre-Omicron (including D614G, Alpha, Beta, Delta, and other non-Omicron variants in South Africa) infection or a primary schedule of vaccination.

• Om denotes immune histories consisting in a single exposure to SARS-CoV-2 via an Omicron infection.

• pOm − pOm denotes immune histories including at least two antigen exposures with first exposure being either pre-Omicron strain or vaccination and second exposure either pre-Omicron strain or vaccination.

• pOm − Om denotes immune history including at least two antigen exposures with the first exposure being either pre-Omicron strain or vaccination and the second exposure being an Omicron infection.

Where S ∈ {pOm, Om} and:

• pOm denotes protection against pre-Omicron strains.

• Om denotes protection against Omicron.

In this study, since the antigenic difference between Omicron and all previously circulating strains (D614G, Alpha, Beta, Delta) is much larger than antigenic differences among previously circulating strains (17, 25, 31), we do not differentiate immunity between pre-Omicron variants (all pre-Omicron variants confer the same type of homologous and heterologous immunity to each other, and we also consider the vaccine to be antigenically similar to pre-Omicron strains since it is based on wild type strain). For simplicity, considering the low vaccination rate in South Africa, we assume that a full vaccination schedule confers similar levels of protection as an infection with a pre-Omicron strain (i.e., D614G, Alpha, Beta, Delta). In reality, a full schedule of vaccination likely confers lower level of protection against infection over long time scales due to waning, irrespective of the vaccine platforms (15). Differences in long-term protection against onward transmission and hospitalization remain unclear. For simplicity and given the relatively short time scale being considered, we also assume that the first two antigen exposures dominate the acquisition of immunity and subsequent infections by Omicron beyond the first two may not change immune memory against Omicron (26). An exhaustive list of all immunity protection scenarios considered in this study is shown in Table S3.

5.4 Estimating the growth advantage of Omicron over Delta during its initial emergence

When Omicron was discovered in South Africa, the Delta epidemic had already declined and the Delta variant was circulating at low level in most locations (20). In Figure S7, the dots show the logarithmic of SARS-CoV-2 weekly incidence in the District of Kenneth Kaunda between weeks 35 and 48 in 2021. The epidemic curve can be viewed as the supposition of an exponential decay and an exponential growth, with the transition occurring around week 45 of 2021 (coinciding with the emergence of Omicron in South Africa). Here we assume that the exponential decay (prior to week 45, 2021) was driven by the Delta variant and the exponential growth (post week 45, 2021) was driven by the Omicron variant. We can thus model the epidemic curve as supposition of exponential decay and exponential growth, i.e.

Here α is the growth rate of Delta and α is the growth rate of Omicron and A andB are the initial Incidence rate for Delta and Omicron when t = 0. We fit this function to the observed epidemic curve between week 35 and week 48 using maximum likelihood method. We find a growth rate of -0.063 per day for the Delta variant (exponential decay) and 0.275 per day for the Omicron variant (exponential growth), indicating a growth advantage of 0.338 per day of Omicron over Delta in Dr Kenneth Kaunda District. Figure S7 shows the results of the fitting.

5.5 State-space transmission model for Delta variant and projection of Delta spread from weeks 35 to 45, 2021

Here we consider a “Susceptible-Infectious-Recovered-Susceptible” (SIRS) model for Delta that that tracks infection history up to 3 repeat infections/immunizations. The equations governing the model are as follows:

With the following expression for the force-of-infection:

The infectious periods for and i₂ can be expressed as follows with respect to i₀’s:

The definitions of the state variables are presented in Table S4, and the definitions of the model parameters are in Table S5. We initialized the model based on the reconstructed exposure history (Figure S6) by week 34. The only free parameter of the model was the rescaling factor on transmissibility α_scaling (Table S5). We optimized α_scaling so that the projected curve of new incidence would decay at a rate of -0.063 day^-1, matching the observation (Figure S7). The best fit for α_scaling was 0.44. We further projected the SARS-CoV-2 antigen exposure history from week 35 to week 45 (prior to the emergence of the Omicron variant (Figure S8).

5.6 State-space transmission model for the Omicron variant and projection of Omicron spread from week 45, 2021 to the end of the Omicron wave

For simplicity, we consider a hypothetical scenario where Omicron has successfully displaced all other circulating variants in South Africa and explore how the transmission dynamics of Omicron is shaped by the immune history of previously circulating strains and vaccination. Accordingly, at the time of writing, Omicron had replaced Delta in many countries that report variant-specific prevalence estimates, including South Africa. We do not consider variant co-circulation and the emergence of new variant during the Omicron wave, although such scenarios are certainly possible. Similar to the Delta variant, we consider a “Susceptible-Infectious-Recovered-Susceptible” model for Omicron that tracks infection history up to 3 repeat infections/immunizations, with additional Omicron-specific properties of immune evasion and enhanced transmissibility. The equations governing the model are as follows:

With the following expression for the force-of-infection:

The infectious periods for and can be expressed as follows with respect to i₀ ‘s:

The definitions of the state variables are presented in Table S6, and the definitions of the model parameters are in Table S7.

Based on the transmission model, we explored how the degree of Omicron’s immune evasion against infection and onward transmission would shape the trajectory of the epidemic (See Section 5.2 and Table S3 for the definition of and We scanned through values of and ranging from 0 to 1, with a step size of 1/30. We also considered a potential change in the generation interval (GI) of Omicron when compared to the Delta variant. We explored possible values for Omicron’s generation time, including 3, 4, 5, and 6 days. For each pair value pair of and and GI, we fit the ratio of basic reproduction number between Omicron and Delta so that the growth rate of Omicron matched the observed initial growth of the Omicron wave (Figure S7). We fit the fraction of individual in the i₀ compartment at week 45, 2021 so that the peak of the projected incidence of Omicron infections matched the observed Omicron case incidence.

We then calculated the characteristics of the projected Omicron wave, including the estimated 1) infection attack rate, 2) epidemic duration, 3) fraction of reinfections/breakthrough infections among all infections, 4) the relative reduction of realized GI (average GI over both primary infections and reinfections/breakthrough of infections) with respect to intrinsic GI, and the 5) infection case ratio (number of cases reported to the Dr Kenneth Kaunda District during the Omicron wave divided by the total number of projected Omicron infections), for a given and and GI.

Figure S9 visualizes projected characteristics of the Omicron wave as a function of and and GI. Figure 4A-E represents a scenario where the GI of Omicron is 4 days, shorter than Delta (11, 17, 18, 26), where white dots represent our best knowledge of the degree of Omicron’s evasion of prior immunity against infection and onward transmission (64).

6.1 Omicron transmission model

For the Omicron transmission model, we first consider the same parameters as in the reference scenario as shown in Figure 4F, with generation time of 4 days, , , and the estimated . However, for the two population groups of interest 1) those who were primed by Omicron 2) those who have experienced a Omicron reinfection/breakthrough we consider their relative risk (with respect to naïve population) of acquiring Omicron infection as and, respectively. In Figure S10A, we use the transmission model to evaluate the growth rate of Omicron when Omicron is reintroduced into the population for all combinations of and ranging from 0 to 1. A growth rate larger than 0 indicates that the Omicron variant is above the epidemic threshold, leading to a recurring fifth epidemic wave, while a growth rate lower than 0 indicates that the Omicron variant will not trigger another outbreak after the fourth wave. We further consider a scenario (Figure S10B) where the contact rate is twice that of the one during the fourth wave, i.e., α_scaling= 0.44 × 2 = 0.88.

6.2 Delta transmission model

To evaluate the risk of Delta recurrence, we consider the same Delta transmission model described in Materials and Methods Section 5.5. However, for the population groups 1) who were primed by Omicron 2) who have experienced a Omicron reinfection/vaccine breakthrough, we consider their relative risk (with respect to naïve individuals) of acquiring Delta infection as and , respectively. In Figure S10C, we use the transmission model to evaluate the growth rate of the Delta variant when it is reintroduced into the population for any combination of and ranging from 0 to 1. We further consider a scenario (Figure S10D) where the contact rate is twice that of the one during the fourth wave, i.e., α_scaling = 0.44 × 2 = 0.88 (Table S5).

6.3 Transmission model of the hypothetical variant X

For the transmission model of a hypothetical new variant X (Figure S10E), we consider that variant X has the same basic reproduction number and generation time as the Delta variant, and the contact rate is twice that of the fourth wave, i.e., α_scaling = 0.44 × 2 = 0.88. We additionally assume that variant X is antigenically equally distinct from both Omicron and pre-Omicron variant so that the relative risks of reinfection are equal irrespective of the primed strain i.e.,. The rest of the model is the same as for Delta.

• Download figure
• Open in new tab

Fig. S1:

Vaccination rate by vaccine types, including J&J/Janssen Ad26.COV2.S (J&J) and the Pfizer/BioNTech BNT162b2 (Pfizer). The dashed lines indicate the end of PHIRT-C (August 28, 2021).

• Download figure
• Open in new tab

Fig. S2:

(A) Weekly incidence (blue bar) and cumulative incidence (red line) of SARS-CoV-2 infection among participants in the rural cohort. (B) Same as (A) but in the urban cohort.

• Download figure
• Open in new tab

Fig. S3:

Distribution of the size of the household infection cluster at different household sizes across 222 households in both the rural and urban cohort. A household infection cluster with size larger than one is defined as a group of infections within the same household with at least two infection episodes within the same cluster with infection time separated no more than 14 days. An isolated infection episode is considered an infection cluster with cluster size one. The sizes of the dots are proportional to the frequency of occurrence for each “household size – cluster size” pair among 192 infection clusters.

• Download figure
• Open in new tab

Fig. S4:

(A) Discrete time survival analysis using piecewise exponential model (Poisson regression). Hazard ratios (HR) along with 95%CIs are reported as solid dots and horizontal lines. The hollow dots are reference class for each of the categorical variable. Comparing to results presented in Figure 3C, we consider a censoring time window from the time of infection to 15 days after viral RNA clearance for each infection episode. (B) Discrete time survival analysis using piecewise logistic model (logistic regression). Comparing to results presented in Figure 3C, we consider the piecewise logistic model rather than Poisson model, with the same censoring time window. Odds ratios (OR) along with 95%CIs are reported as solid dots and horizontal lines. The hollow dots are reference class for each of the categorical variable. *Indicates p<0.05; ** indicates p<0.01; *** indicates p < 0.001. Abbreviations: HIV- (HIV-uninfected individuals), PLWH+ CD4 <200 (Persons living with HIV, CD4+ T cell count under 200 cells/ml), PLWH+ CD4 >=200 (Persons living with HIV, CD4+ T cell count equal or above 200 cells/ml).

• Download figure
• Open in new tab

Fig. S5:

(A) Same as Figure 4F but for a low immune escape (LE) scenario with and . (B) Same as Figure 4F but for a high immune escape (HE) scenario with and .

• Download figure
• Open in new tab

Fig. S6:

SARS-CoV-2 antigen exposure history by the end of PHIRST-C (September 2021) at the District of Dr Kenneth Kaunda.

• Download figure
• Open in new tab

Fig. S7: Estimated case growth rate of Delta and Omicron between weeks 35 and 52 of 2021.

Dots are the logarithmic of the weekly incidence of SARS-CoV-2 cases reported to the District of Dr Kenneth Kaunda between weeks 35 and 52 of 2021. The blue line is the fitted exponential decaying of Delta Incidence while the red line is the fitted exponential growth of Omicron wave and the dashed line is the fitted convolution of the two variants.

• Download figure
• Open in new tab

Fig. S8:

Projection of SARS-CoV-2 antigen exposure history from the end of PHIRST-C (week 35 of 2021, dashed line) and until the emergence of Omicron (week 45 of 2021, Figure S7) at the District of Dr Kenneth Kaunda.

• Download figure
• Open in new tab

Fig. S9: The characteristics of the Omicron epidemic wave as a function of and and generation time.

Each row corresponds to different values of generation time (GI) of primary infection. Column one corresponds to ratio of basic reproduction number between Omicron and Delta. Column two corresponds to infection attack rate. Column three corresponds to duration of the epidemic. Column four corresponds to fraction of reinfections/breakthrough of infections among all infections. Column five corresponds to the relative reduction of realized GI (average GI over both primary infections and reinfections/breakthrough infections) with respect to the intrinsic GI.

• Download figure
• Open in new tab

Fig. S10:

Possible post-Omicron futures, exploring potential resurgences of Omicron, Delta, or a new variant X, under different contact scenarios. Projections are based on the reconstructed immune histories at the end of January 2022 shown in Fig 4F (i.e., given a reference scenario for Omicron’s immune escape and ). (A) Risk of recurrence of Omicron: phase diagram of the growth rate of Omicron in a recurrent wave, as a function of the level of protection conferred by Omicron primary infections against Omicron , and the level of protection conferred by Omicron reinfections/breakthroughs against Omicron . Contact rates are assumed to remain the same as during the Omicron wave. (B) same as (A) but assuming the contact rate is twice of that during the Omicron wave. (C) Risk of recurrence of Delta: Phase diagram of the growth rate of Delta, after the initial Omicron wave has subsided, as a function of the level of protection conferred by Omicron primary infection against Delta , and the level of protection conferred by Omicron reinfections/breakthroughs against Delta . Contact rates are assumed to remain the same as during the Omicron wave. (D) same as (C) but assuming the contact rate is twice of that during the Omicron wave. (K) Risk of occurrence of hypothetical new variant X, where X is at equal antigenic distance of Delta and Omicron: Phase diagram of the growth rate of variant X, after the initial Omicron wave has subsided, as a function of the level of protection conferred by any variant primary infection on infection with , assuming Omicron and pre-Omicron infections confer the same level of protection against X), and the level of protection conferred by Omicron reinfections/breakthroughs on X infection,. Contact rate is assumed to be twice of that of the Omicron wave.