Estimating the time-dependent effective reproduction number and vaccination rate for COVID-19 in the USA and India

1.   Introduction

A disease outbreak (COVID-19) of an uncommon pneumonia triggered by a new coronavirus, SARS-CoV-2, appeared in Wuhan, China, in December 2019 ^[1], has instigated at least 631, 989, 826 confirmed cases, with 6, 578, 767 deaths in nearly 230 nations and territories by October 21, 2022 ^[2]. The COVID-19 crisis was acknowledged as a pandemic by the World Health Organization (WHO) on March 11, 2020. The earliest case and death in India were reported on January 30, 2020 and March 13, 2020, respectively ^[3]. In the early outbreak, pharmaceutical interventions, for instance, drugs and vaccination, were unavailable, thereafter controlling the pandemic greatly depends only on non-pharmaceutical interventions, such as lockdown, maintaining social distancing, contact tracing, testing, using face masks, etc ^[4,5,6,7]. However, in India, the vaccination campaign with a single and double dose was started on January 15, 2021 and February 13, 2021, respectively. The double-dose vaccination campaign started on December 13, 2020 in the USA ^[8].

The estimation of the reproduction number of a disease is particularly important in modeling the transmission of infection as well as notifying and assessing control policies ^[9,10]. On the other hand, the instantaneous (or effective) reproduction number $R_t$ exposes the limit of disease spread with the impact of an intervention or population immunity. Public health officials and policymakers use $R_t$ to evaluate the efficacy of control measures and several other factors that influenced the disease transmission at certain points of time and to update policy accordingly. $R_t$ estimates could be exercised to evaluate the effectiveness of control measures or as a near real-time indicator of epidemic progress. For these intents, estimates need to be precise and accurately signify uncertainties, and for near real-time scrutinizing, they also need to be timely. Consequently, $R_t$ estimation is crucial for assessing the achievement of public measures. However, $R_t$ estimation from available data offers numerous challenges, with crucial inferences for the interpretation of the course of the epidemic ^[11,12]. Though, the parameter assumptions and model structure makes the estimation of $R_t$ sensitive ^[11]. Poor quality or unavailability of data, on the other hand, often impedes the utilization of several estimation methods, for instance, serial intervals of a disease that are typically required for the estimation of $R_t$ ^{[12,13,14,15]}.

Specific hypotheses and data estimates are required for computing the exact value of $R_t$. Several methods are used to estimate the effective reproduction number and time-dependent rate in the literature. Cazelles et al. ^[16] utilized stochastic models for the transmission dynamics of a disease coupled with particle Markov Chain Monte Carlo Algorithm and estimated the time-dependent transmission rate. Nishiura et al. ^[17] discussed a likelihood-based method for the estimation of $R_t$ from the data of early epidemic growth. Bettencourt and Ribeiro ^[18] estimated $R_0$ and $R_t$ by using the compartmental Susceptible-Infectious-Recovered (SIR) model and incidence data of H5N1 influenza in humans. Pollicott et al. ^[19] and Kong et al. ^[20] introduced an inverse method for epidemic models that extract the time-dependent transmission rate from the incidence or prevalence data. Wang et al. ^[21,22] developed an inverse method and estimated the time-dependent transmission rate and death rate using a generalized boosting model (GBM) and a mechanistic ordinary differential equation (ODE) model. Hasan et al. ^[23] developed a novel method to estimate the effective reproduction number based on the Kalman filter and a low pass filter method. The authors used a simple SIRD (Susceptible-Infectious-Recovered-Dead) model and augmented the model into a discrete-time stochastic model. This method is quite beneficial since it does not require serial intervals, which forms the process easier without falling the estimation quality.

In this paper, we estimate the time-dependent effective reproduction number and vaccination rate using the method developed by Hasan et al. ^[23] for COVID-19 in the USA and India. We use an augmented SVEIR model and the data of daily new cases, cumulative cases, daily new deaths, cumulative deaths, and fully vaccinated individuals concerning COVID-19.

2.   A discrete-time stochastic augmented SVEIR model

We estimate the effective reproduction number, $R_t$, and time-dependent vaccination rate, denoted by $\xi_t$, by using a simple compartmental SVEIR model. Figure 1 shows the scheme of the model. Let us denote $S, V, E, I,$ and $R$ respectively as susceptible, vaccinated, exposed (not yet infectious), active/confirmed, and recovered compartment. Let $N = S+V+E+I+R$, where $N$ denotes the total population.

In the model (2.1), the susceptible individuals are vaccinated at rate $\xi_t$, which is time-dependent. $\beta_t$ signifies the average number of contacts per individual per unit time and is time-dependent due to interventions. Since the vaccine is imperfect, vaccinated individuals can also get the infection, but at a reduced rate, $\sigma \beta_t$, where $1-\sigma$ characterizes the vaccine efficacy. The immunity acquired by vaccination wanes with time, and vaccinated individuals can become susceptible again at a rate $\phi$. $a$ stands for the transition rate of exposed individuals becoming infectious, $\alpha$ denotes the recovery rate of infectious individuals, and $\delta$ is the disease-induced death rate. $\gamma$ represents the antibody waning rate since the immunity acquired by infection also wanes with time. Hence, the model is written in the following non-linear ordinary differential equations:

To use the method ^[23], we also need some more information about the Case Fatality Rate (CFR) and the average infectious time, $T_i$, and

For COVID-19, we consider $T_i = 14$ ^[21], however it may vary from 1.5 to 30 days amongst distinct diseased persons ^[24]. CFR is assumed to be 1%.

Covaxin and Covishield, two COVID-19 vaccines were advanced and being utilized in India. Bhatnagar et al. ^[25] studied that the vaccine efficacy of full (two doses) vaccination was 71% (95% CI: 57-81%) with Covaxin/BBV152 and with covishield/AZD1222, it was 85% (95% CI: 79-89%). The vaccine efficacy was highest at 94% (95% CI: 86-97%) for 6-8 weeks between two doses of Covishield and 93% (95% CI: 34-99%) of Covaxin. In the USA, the COVID-19 mRNA vaccines BNT162b2 (Pfizer-BioNTech) and mRNA-1273 (Moderna) COVID-19 vaccines have been completely permitted by the Food and Drug Administration (FDA) for individuals aged 18 years and older and individuals aged 16 years and older, respectively ^[26]. In clinical trials, mRNA based COVID-19 vaccines, BNT162b2 (Pfizer Inc, New York, USA; BioNTech Manufacturing GmbH, Mainz, Germany; Comirnaty) and mRNA-1273 (Moderna Inc, Cambridge, USA; Spikevax) were greatly effectual (95% and 94%, respectively) against symptomatic COVID-19 infection ^[27]. In our estimation, vaccine efficacy is to be assumed 95% and 85% for the USA and India, respectively.

To compute the effective reproduction number, $R_t$, we can still apply the next generation matrix method ^[28,29]. According to Lemma 1 in ^[30], we define

Following ^[28,29], the effective reproduction number is defined as the spectral radius of the next generation matrix ($FV^{-1}$), i.e.,

The approximation is under the assumption that government intervention is taken at an early stage so that $S(t)+\sigma V(t) \approx N$. This is the case especially for the emerging diseases.

We further augment the SVEIR model (2.1) by considering the following seven equations:

where the first equation takes into account the cumulative cases ($C$), second equation provides the daily number of new cases ($D$), third equation gives cumulative deaths ($CD$), fourth equation offers the daily number of new deaths ($DD$), fifth equation produces the cumulative number of vaccinated individuals ($P$), the sixth and seventh equations say that the effective reproduction number ($R_t$) and vaccination rate ($\xi_t$) are assumed to be piecewise constant functions with one day time interval jump.

Further, by discretization and applying the forward Euler method, the discrete-time augmented SVEIR model is given as follows:

We estimate the basic reproduction number, $R_t$, by using the cumulative cases, new daily cases, cumulative deaths and new daily deaths, and vaccination rate, $\xi_t$, by using the full (with two doses) vaccination data. The data has been collected from the World Health Organization (WHO) website https://covid19.who.int/ and https://ourworldindata.org/coronavirus on a daily basis. Since the frequency of the data could be once a day which is low, then by utilizing a modified Akima cubic Hermite interpolation, the new reported data could be interpolated and it fits with $\Delta t$ time steps. The $\Delta t$ time step is taken as $0.01$ in our simulation, i.e., 100 times discretization within one day time interval. For simplicity, we outline the augmented state vector as follows:

and the discrete time augmented SVEIR model (2.11) can be rewritten in the following form:

where $g$ is the nonlinear expression of the equations in the right-hand side of the augmented SVEIR model (2.11) and $v$ denotes the uncertainty to the inaccuracies in the model owing to interpretation in the modeling. The uncertainty is supposed to be a zero mean Gaussian white noise with known covariance $Q_F$. It is to make the computation simple since the real epidemic data generally follow Gamma distribution. $Q_F$, in practice, could be taken as a tuning parameter for the EKF. Hence, the model becomes in the form of a discrete-time stochastic augmented SVEIR model.

New daily cases, cumulative cases, new daily deaths, cumulative deaths, and cumulative vaccinated individuals can be incorporated into the model including the ensuing output vector

where $z$ signifies uncertainties owing to the incorrect testing outcomes. It is also supposed the uncertainties to be a zero mean Gaussian white noise with known covariance $R_F$. Similar to $Q_F$, $R_F$ could also be taken as a tuning parameter. Ensuing the convenient data, the measurement/data matrix $C$ is considered to be

A two-stage filtering approach: We use a two-stage filtering approach for the estimation of the time-dependent effective reproduction number, $R_t$, and vaccination rate, $\xi_t$. This method consists a low-pass filter and the EKF. The EKF is an expansion of Kalman filter (KF) for nonlinear systems. The KF is centered on Bayesian estimation and is an optimal linear filter. The EKF is established based on the linearization of the non-linear system around its estimate. Because of this linearization, the stability and optimality of the EKF cannot be ensured. However, the EKF may provide a fairly good estimation if the non-linearity is not very extreme. In the estimation process, the state variables, $R_t$ and $\xi_t$ are estimated by using the EKF under uncertainties in the number of new reported cases which could be a consequence of the delays in reporting.

Let $\hat{y}(l)$ be a vector state, estimated by the EKF. Operating the first-order Taylor series expansion to $g$ at $\hat{y}(l)$, which yields

Here $J_g (\hat{y}(l))$ denotes the Jacobian matrix of $g$, specified by

with

The EKF contains two stages: forecast and update. The model (2.11) is utilized to forecast the upcoming state and covariance and upgrade them after attaining measurement/new data. The EKF could be pondered as one of the straightforward dynamic Bayesian networks. EKF recursively computes the estimate of correct state values over time utilizing a mathematical process model and next measurements while recursive Bayesian estimation computes estimates of an unidentified probability density function over time utilizing a mathematical process model and incoming new data ^[31]. Let $\hat{y}(m|n)$ signifies the estimation of $y$ at time $m$ considering observations up to and together with the time $n \leq m$. The Kalman filter algorithm is given as follows.

Forecast:

Update:

where $p(l|l)$ signifies a posteriori estimate covariance matrix. In the step, "update", a low pass filter approach concerning a rational transfer function is utilized to detach short term fluctuation at time step $l$, and is specified by

where $x_n$ is a window length along the data. We take $x_n = 3/\Delta t$ for our case.

To ensure the estimation quality, we compute a Relative Root Mean Square Error (RRMSE) between the reported and estimated cases. The RRMSE is expressed by

where $O_d$ denotes the total number of observed time in days and $Y \in \left\lbrace C, D, CD, DD, P \right\rbrace$.

We illustrate the observed and estimated numbers of COVID-19 cases in India in Figure 2 and obtained a good agreement between the data and estimated numbers. We also plot the reported and estimated number of fully vaccinated individuals in Figure 3. In the estimation procedure, the observation and process covariance matrices are chosen as tuning parameters and are taken as $R_F = diag(100\; 10\; 5\; 1\; 1)$ and $Q_F = diag(100\; 10\; 10\; 10\; 5\; 5\; 5\; 5\; 5\; 0.5\; 0.5\; 0.5)$, respectively. These tuning parameters are attained from error and trial and adopted such that RRMSE between the reported and estimated data are adequately small. For our case, the RRMSEs are given in Table 2. Here, we notice that the obtained RRMSEs are sufficiently small. The error could be reduced more by choosing distinct values of $R_F$ and $Q_F$.

For the USA, the observed and estimated COVID-19 and vaccinated cases are plotted in Figure 6 and 7, respectively. The observation and process covariance matrices are chosen as $R_F = diag(1000\; 100\; 50\; 10\; 10)$ and $Q_F = diag(100\; 10\; 10\; 10\; 5\; 5\; 5\; 5\; 5\; 0.5\; 0.5\; 0.5)$, respectively. The RRMSEs for the USA are also given in Table 2, which are relatively smaller than India.

We finally plot the time-dependent effective reproduction number $R_t$ and vaccination rate $\xi_t$ for India in Figure 4 and 5, respectively, and for the USA in Figure 8 and 9, respectively. We do not take average out of the new cases data, but consider the likelihood estimation of a new case at one day to depend also on the estimation of the previous three days. The estimated curves of $R_t$ and $\xi_t$ for the USA and India are also showing spikes and serrations with data.

We also provide the short-term forecast for India from August 23, 2022 to December 31, 2022. The forecasting scenario is shown in Figures 2, 3, 4, and 5 in green colored curves. With the current scenario, the prediction shows that the new daily cases and deaths are steadily decreasing in India. However, the effective reproduction number will be greater than one by December 31, 2022.

The predictions for the USA from August 17, 2022 to December 31, 2022 are also provided in Figures 6, 7, 8 and 9 in green colored curves. These results also show that the number of new daily cases and deaths is decreasing. From Figure 8, we observe that the effective reproduction number will be greater than one by December 31, 2022.

3.   Discussion

Several mathematical models and approaches have been evolved and utilized to estimate various types of time-dependent rates and reproduction numbers. Estimating the time-dependent vaccination rate with the effective reproduction number for COVID-19 progress is essential for making public health decisions to regulate the pandemic. We here estimate the time-dependent effective reproduction number $R_t$ and vaccination rate $\xi_t$ based on model (2.11) with respect to the COVID-19 data in the USA and India by using a low pass filter approach and EKF developed by Hasan et al. ^[23]. Our estimated $R_t$ and $\xi_t$ chase the dynamics of the model, showing the spikes and serrations with data. The prediction results show that the new daily cases and deaths are decreasing in the USA and India. We observed that the effective reproduction number would be greater than one by December 31, 2022 in both countries. Which implies that it is still necessary to maintain safety and preventive measures to control COVID-19. Our results are beneficial for the policy/health makers to track the status of the effective reproduction number that whether it is greater or less than one.

The inverse method, given by Pollicott et al. ^[19] and Kong et al. ^[20], needs an explicit formula to estimate the time-dependent rates. In comparison, the method we used in this work is more straightforward and has various advantages. From the modeling perspective, the time-dependent effective reproduction number and vaccination rate follow the dynamics according to the data. This method does not need serial intervals' information. In this case, the provided time series data have drastic changes or instead do not have many changes. The resulting time-dependent rates show the same serrations and spikes. Consequently, the latest information could be elaborated from the data dynamics.

Various rates in epidemiological models for contagious diseases vary significantly over time. It is further feasible to estimate other time-dependent rates, such as immunity waning rate and vaccine efficacy, that decrease over time. The model we utilized in this work could also be applied to COVID-19 transmission dynamics in any other country, specifically when vaccination is employed. We have not considered reporting ratio in our model, however, by considering it, time-varying estimates of COVID-19 infection will be observed in future. Using an extended multi-strain SVEIR model (including multiple infectious compartments with different strains), estimation of the strain-specific effective reproduction number will be addressed in the future.

Acknowledgments

The research work of Sarita Bugalia is supported by the Council of Scientific & Industrial Research (CSIR), India [File No. 09/1131(0025)/2018-EMR-I]. The research work of Jai Prakash Tripathi is supported by the Science and Engineering Research Board (SERB), India [File No. ECR/2017/002786]. The research work of Hao Wang is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). We are highly thankful to all the anonymous reviewers and editor for their insightful comments and suggestions, which helped us to improve the manuscript considerably.

Conflict of interest

The authors declare there is no conflict of interest.