Characterizing the effective reproduction number during the COVID-19 epidemic: Insights from Qatar’s experience

INTRODUCTION

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic is the most serious global health challenge in recent history [1,2]. Coronavirus Disease 2019 (COVID-19) morbidity and mortality has imposed unparalleled burdens on healthcare systems worldwide, and necessitated unprecedented restrictions on mobility and on social and economic activities [3,4]. Tracking and monitoring each wave of infection have become essential to avoid the adverse consequences of infection transmission [5-8]. With such serious consequences to the healthcare system, economy, and society, decisions regarding the escalation or easing of restrictions have become a critical facet of policymaking since the discovery of the virus in December of 2019 [6,9,10].

The effective reproduction number (R_t), the average number of secondary infections each infection is generating at a given point in time [6,11-13], has proven to be one of the most influential tools in monitoring and tracking the epidemic, and informing the escalation and easing of restrictions [6,11-13]. In the present study, we report on the use of R_t in Qatar, to understand epidemic dynamics and to establish national policy decisions and public heath interventions, in what has become a successful country experience.

Qatar is a peninsula in the Arabian Gulf with a diverse population of 2.8 million people [5,14] that has been affected by three SARS-CoV-2 epidemic waves [5,6,15-19]. The first wave started with the introduction of the virus in February of 2020 and peaked in late May 2020 [5,6]. The second wave started in mid-January, 2021, and was triggered by the introduction and expansion of the B.1.1.7 (Alpha [20]) variant [15-19,21]. This wave peaked in the first week of March, but was followed immediately by a third wave that was triggered by introduction and rapid expansion of the B.1.351 (Beta [20]) variant, which started in mid-March and peaked in mid-April, 2021 [15-19,21].

Two forms of R_t estimation have been used in Qatar to inform the national COVID-19 response, and each proved to have its own intrinsic public health value. The first is the real-time “empirical” estimation which is done by calculating R_t directly from diagnosed cases. Different methods were explored for estimating the empirical R_t (henceforth, ), and based on this exploration the Robert Koch Institute method [13,22] was used for feasibility, ease of use, and functionality. The second estimation method was model-based by calculating R_t using a population-level compartmental transmission dynamics model [6,23], hereafter designated as .

METHODS

RESULTS

The calculated using the Robert Koch Institute method captured effectively the evolution of the epidemic through its three waves, starting from the first wave (the wild-type variant wave) [5,6], the second (B.1.1.7) wave [15-19,21], and the third (B.1.135) wave [15-19,21] (Figure 1A). It also captured and correlated with key response landmarks, such as partial lockdowns during the three waves and subsequent easing of public health restrictions, and major social events that led to transient increases in the social contact rate in the population. It further captured transient fluctuations that were associated with significant clusters of infection, especially during low-incidence phases between August 1, 2020 and January 17, 2021, and between May 25, 2021 and August 18, 2021.

• Download figure
• Open in new tab

Figure 1.

Effective reproduction numbers and in Qatar. A) Trend in and , April 1, 2020 to August 18, 2021, and association with major events, response landmarks, and introduction and expansion of the B.1.1.7 and B.1.135 variants. B) Trend in for only the B.1.1.7 variant cases, February 1, 2021 to April 1, 2021. was estimated using the Robert Koch Institute method [22] applied to symptomatic case series data. The dashed green line represents the threshold of R₀ = 1.

The epidemic expansion of B.1.1.7 cases starting from January 18, 2021 was associated with a large and rapid increase in (Figure 1A), suggesting the higher infectiousness of this variant. calculated using only B.1.1.7 case series data is shown in Figure 1B and demonstrated higher values, confirming further the higher infectiousness of this variant. for only the B.1.1.7 variant averaged 1.45 during the exponential growth phase of the second (B.1.1.7) wave (February 1-22, 2021). It was unstable during the first two weeks of this wave (January 18-31, 2021; not shown), as transmission appears to have been influenced by one or more superspreading events that were not representative of the average community transmission. It was also unstable after April 1, 2021, as the number of daily B.1.1.7 cases was too small.

The first sensitivity analysis on the estimated , in which the time series of all diagnosed cases was used instead of only symptomatic cases, showed overall excellent agreement regardless of the input-data source used to calculate (Figure 2A). The Pearson correlation coefficient was 0.914 (p-value<0.001) and the Spearman correlation coefficient was 0.913 (p-value<0.001), both confirming the excellent agreement. There were only few noticeable differences that occurred when the number of diagnosed cases was too small (periods when SARS-CoV-2 incidence was low); thus, was more susceptible to transient variation in the number of diagnosed cases, such as due to sporadic, random PCR testing campaigns. Peaks in were also slightly larger using only symptomatic cases versus all diagnosed cases.

• Download figure
• Open in new tab

Figure 2:

Sensitivity analyses of estimated using the Robert Koch Institute method. A) Sensitivity analysis using the time series of all diagnosed cases instead of only symptomatic cases in estimating . B) Sensitivity analysis using the time series of hospital admissions in acute-care beds instead of symptomatic cases in estimating . C) Sensitivity analysis using the time series of hospital admissions in ICU-care beds instead of symptomatic cases in estimating . D) Sensitivity analysis using different values for the generation time in estimating . The dashed green line represents the threshold of R₀ = 1.

The second sensitivity analysis, in which the time series of acute-care hospital admissions was used to proxy the epidemic trend, instead of the time series of symptomatic cases, showed rather strong correlation between estimates (Figure 2B). The Pearson correlation coefficient was 0.512 (p-value<0.001) and the Spearman correlation coefficient was 0.716 (p-value<0.001), both confirming strong agreement. The third sensitivity analysis, in which the time series of ICU-care hospital admissions was used to proxy the epidemic trend, instead of the time series of symptomatic cases, also showed rather strong correlation between the estimates, but also large fluctuations in (Figure 2C). The Pearson correlation coefficient was 0.589 (p-value<0.001) and the Spearman correlation coefficient was 0.550 (p-value<0.001), both confirming strong agreement, but inferior to that for acute-care hospital admissions (Figure 2B versus Figure 2C).

The fourth sensitivity analysis, in which different values for the generation time τ_G were used, showed also strong agreement between different estimates (Figure 2D). The Pearson correlation coefficient was 0.901 (p-value<0.001) and the Spearman correlation coefficient was 0.900 (p-value<0.001), both confirming excellent agreement. The main differences between the estimates occurred in the timing and magnitude of peaks of the epidemic waves, as expected, since variation in generation time changes the rate of epidemic growth [51].

estimated using the Robert Koch Institute method (Figure 3A), Cislaghi method (Figure 3B), Systrom-Bettencourt and Ribeiro method (Figure 3C), Wallinga and Teunis method (Figure 3D), and Cori et al. method (Figure 3E), all showed similar results and were able to capture the evolution of epidemic waves and transient variations due to national public-health response landmarks and major social events. There were also overall strong correlations between them (Table 1 and Figure 4). However, the Systrom-Bettencourt and Ribeiro method (Figure 3C) tended to provide something of an average and was not as sensitive to transient changes in (Figures 3 and 4).

• Download figure
• Open in new tab

Figure 3:

Trend in in Qatar, April 1, 2020 to August 18, 2021, using the A) Robert Koch Institute method [22], B) Cislaghi method [34], C) Systrom-Bettencourt and Ribeiro method [12,37-39], D) Wallinga and Teunis method [35], and E) Cori et al. method [31]. The figure includes the 95% uncertainty or credible interval, as applicable for each method. The dashed green line represents the threshold of R₀ = 1.

• Download figure
• Open in new tab

Figure 4:

Pairwise comparison between estimated using the Robert Koch Institute method [22] and that estimated using the A) Cislaghi method [34], B) Systrom-Bettencourt and Ribeiro method [12,37-39], C) Wallinga and Teunis method [35], and D) Cori et al. method [31]. The dashed green line represents the threshold of R₀ = 1.

There were differences in how rapidly each method detected a change in epidemic dynamics (Figure 4). The Wallinga and Teunis method was the fastest at detecting a change, while the Robert Koch Institute method was the slowest, leading to weaker Pearson and Spearman correlation coefficients between them (Table 1). For instance, the surge in during the first Eid al-Adha after epidemic onset (a festival that occurred between July 30 and August 6, 2020 and is associated with celebrations and social gatherings) was detected on August 1, August 7, August 8, August 11, and August 13 using the Wallinga and Teunis, Cislaghi, Systrom-Bettencourt and Ribeiro, Cori et al., and Robert Koch Institute methods, respectively.

Uncertainty intervals around the estimates of the various methods were narrow overall, except when the number of diagnosed symptomatic cases or the number of PCR tests was small, specifically during the low-incidence phases of the epidemic (Figure 3). Overall, the uncertainty in estimates did not impact the interpretation of the results (Figure 3). The only exception was for the Systrom-Bettencourt and Ribeiro method, as it showed rather wide uncertainty intervals compared to the point estimates for (Figure 3C).

The correlated strongly with the using different methods (Table 1), and provided somewhat of an average of the (Figure 1). For example, and averaged 1.15 and 1.14 during the first wave, respectively. While it captured the three epidemic waves, it could not capture the transient fluctuations in , nor the effects of significant clusters during low-incidence phases.

DISCUSSION

Two forms of R_t estimation were used in Qatar to inform the national COVID-19 response, the “empirical” estimation, , and the model-based estimation, . Both forms proved useful in real-time tracking of epidemic trends, understanding epidemic dynamics, and setting interventions to control transmission, such as application or easing of public health restrictions. Both forms were integral to the national public health response and to formulating evidence-based policy decisions to minimize the epidemic’s toll on health, society, and the economy throughout the phases of this epidemic.

effectively captured the evolution of the epidemic during its three waves, the effects of the response landmarks, such as the partial lockdowns and easing of public health restrictions, and the major social events that affected the social contact rate in the population. Even transient fluctuations in infection transmission that occurred because of significant infection clusters were captured by . Strikingly, the introduction and expansion of the B.1.1.7 variant [21], that resulted in the second epidemic wave, was discovered immediately through monitoring, as there was a sudden large, sustained increase in R_t that coincided precisely with a rapidly increasing number of S-gene “target failures” in PCR testing, even before viral genome sequencing was conducted to confirm the presence and expansion of this variant in the population [18].

While provided an average R_t that closely tracked the average , it did not have the resolution to capture transient changes in , other than major changes associated with the three epidemic waves. Still, was useful and influential, as it was, along with the model that generated it [6,23], the basis for forecasting and future planning, such as forecasting the epidemic time-course and epidemic potential, forecasting healthcare needs of acute-care and ICU-care bed hospitalizations, predicting the impact of social and physical distancing restrictions, planning for easing of restrictions, and forecasting the impact of different mass vaccination strategies [6,23]. Therefore, both forms of R_t complement each other and should be part of any effective COVID-19 national response.

estimation proved robust in sensitivity analyses conducted to assess its utility. Baseline estimation of was based on the time series of symptomatic cases as a proxy of the actual incidence of SARS-CoV-2 infection in the population, which is unknown. Using the time series of all diagnosed cases instead of just symptomatic cases did not appreciably impact estimation, even though PCR testing volume and strategies varied throughout the epidemic. Using the time series of acute-care hospital admissions instead of the time series of symptomatic cases also led to comparable estimates for . This was also the case, but with weaker agreement, when the time series of ICU-care hospital admissions was used to proxy trends in infection incidence. This is not surprising as there is a long delay between onset of infection and ICU-care hospital admission, and the number of ICU-care admissions was relatively small with the low COVID-19 severity in Qatar’s predominantly young and working-age population [5,48]. Variations in the assumed value for the generation time in the estimation did not heavily impact estimates. These findings support the robustness of the approach employed to estimate .

Examination of different methods to estimate demonstrated consistency of the results, generally strong correlations between the estimates, and an acceptable level of uncertainty in them. The only exception was the Systrom-Bettencourt and Ribeiro method which tended to provide something of an average . It was not as sensitive to transient changes in , and had wide uncertainty intervals compared to the point estimates. There were also differences in how rapidly each method detected a change in epidemic dynamics. The Wallinga and Teunis method was the fastest to detect a change, while the Robert Koch Institute method was the slowest. Yet, overall, these findings support the robustness of using these methods in estimation and to guide COVID-19 national responses.

This study has limitations. The estimated and are contingent on the validity and generalizability of input data. The uncertainty/credible intervals estimated here accounted for the uncertainty arising from sampling variation, or from our imperfect knowledge of specific epidemiological quantities, such as the serial interval, but did not account for other sources of uncertainty, such as our imperfect knowledge of the true incidence of infection in the population. To reduce bias due to variation in volume and strategies of PCR testing over time, was calculated using the time series of symptomatic cases, but the distribution of the delay between onset of infection and onset of symptoms may bias these estimates. was estimated using a deterministic compartmental model, but this type of model may not be representative of stochastic transmission dynamics, particularly when the number of infections is small. Despite these limitations, and were able to capture the evolution of the epidemic through its several waves, and to effectively inform the national response and policy decision-making.

CONCLUSIONS

R_t estimations played a critical and influential role in the COVID-19 national response in Qatar. Even though surveillance data of SARS-CoV-2 infection are imperfect and prone to bias, R_t estimations were robust and generated consistent results regardless of the data source used, or the method employed in generating estimates. These findings affirm the value and complementarity of using both and to track the epidemic in real-time and to inform public health decision making at a national level. This can also be done despite low-resource demands, as R_t estimation utilizes existing surveillance data. Moreover, application of some of the estimation methods is feasible even without established expertise in infectious disease modeling. Since the choice of estimation method does not impact the estimates, each country may decide on the best approach, method, and source of data to be used in the estimation, weighing feasibility, ease of use, and functionality, given its specific circumstances.

Author contributions

RB and HHA constructed, coded, and parameterized the mathematical models, and conducted the analyses. HC conducted statistical analyses. HHA and LJA. conceived and led the design of the study, construct and parameterization of the mathematical models, and wrote the first draft of the article. All authors contributed to discussion and interpretation of the results and to the writing of the manuscript. All authors have read and approved the final manuscript.

Supplementary Material

• Download figure
• Open in new tab

Figure S1.

A conceptual diagram illustrating the basic structure of the deterministic mathematical model developed to describe SARS-CoV-2 transmission dynamics in the population of Qatar. The detailed structure of this model and its description are found in [1,2]. In this figure, solid lines denote progression or forward movement from one population compartment to the next, while dashed lines denote backward movement from the present population compartment to the previous one.