Bayesian Framework for Multi-Wave COVID-19 Epidemic Analysis Using Empirical Vaccination Data

Abstract

The COVID-19 pandemic has highlighted the necessity of advanced modeling inference using the limited data of daily cases. Tracking a long-term epidemic trajectory requires explanatory modeling with more complexities than the one with short-time forecasts, especially for the highly vaccinated scenario in the latest phase. With this work, we propose a novel modeling framework that combines an epidemiological model with Bayesian inference to perform an explanatory analysis on the spreading of COVID-19 in Israel. The Bayesian inference is implemented on a modified SEIR compartmental model supplemented by real-time vaccination data and piecewise transmission and infectious rates determined by change points. We illustrate the fitted multi-wave trajectory in Israel with the checkpoints of major changes in publicly announced interventions or critical social events. The result of our modeling framework partly reflects the impact of different stages of mitigation strategies as well as the vaccination effectiveness, and provides forecasts of near future scenarios.

1. Introduction

Throughout a typical outbreak of COVID-19 in a certain country, the first stage of the pandemic is most likely to have a trajectory with exponential growth [1], which is plausibly explained by commonly used statistical models, e.g., the epidemic growth models, such as the Gompertz model [2], or compartmental models [3,4].

However, due to the fact that the COVID-19 pandemic has been hitting most countries worldwide with various patterns of epidemic development for over one year, multiple waves of infection growth have emerged after the initial phase; see [5]. Observing a complete epidemic trajectory of COVID-19 in a certain country, what we could commonly expect is an initial wave of exponential growth followed by multiple waves seperated by corresponding change points. The change points indicate the nationwide changes of effective factors, including mitigation measures [6], containment policies, groups of imported cases, social movement and large-scale clustering events, increasingly infectious variants for pandemics, etc. Considering the complexities, a typical deterministic model with a fixed set of parameters could hardly handle the explanatory analysis of the complete epidemic trajectory due to the limitation of real data and the shifting dynamics corresponding to change points. Hence, the integrated modeling class, such as piecewise models, will unfold the functionalities adapted to trajectories featuring long-term study periods.

Vaccination is another essential factor that considerably affects the transmission of COVID-19, especially for the recent circumstances in numerous countries which have accomplished millions of injections of vaccine doses; see [7]. Unlike most of the other pandemics in the past, thanks to the open source of empirical vaccination data, the exact numbers of daily vaccinations can be directly implemented in the model instead of being approximated through the complex dynamics. Though the vaccination has not produced any impact until the inoculations facilitated by the government, which normally began at the latest phase of the epidemic [8], it enables a substantial decline on susceptible individuals; see [9,10]. Thanks to the high protection rate [11] and the continuously increasing population of vaccinated individuals, the vaccination has brought an obvious influence on the long-term study of the epidemic trajectory.

The ODE-based compartmental models, such as the susceptible–infected–recovered (SIR) model [12,13], play a dominant role in the studies of epidemic outbreaks for infectious diseases. This typical class of methodologies, which simulates population compartments and their shifting ratios from one to another (e.g., the susceptible individuals become infectious), mainly focuses on the estimation and prediction of the epidemic trajectories and effectiveness of containment measures. Based on the idea of depicting the multi-wave trajectory of a country with a highly vaccinated population, in our study, we establish a parsimonious dynamic system of the susceptible–exposed–infected–recovered (SEIR) model [14,15,16] which is modified by piecewise modeling and the additional usage of empirical vaccination data. Parameter inference that competently reflects the corresponding trajectories remains a critical issue for the model calibration. Among them, the transmission rate is particularly sensitive to epidemic changes emerged at different checkpoints throughout the full study period. Hence, the Bayesian framework [17] is utilized for the inference of parameters, due to the complexities and insufficient data of modeling.

Various Bayesian methods have been employed for the epidemiological and related models for studying COVID-19, including the Bayesian MCMC algorithm for the inference of statistical and epidemiological parameters [18,19,20,21], Naïve Bayes method for feature classification in the diagnosis of COVID-19 [22], Bayesian network for severity classification [23,24], spatial geo-referencing analysis [25,26], etc. A few studies have already presented remarkable attempts on methodologies that integrate Bayesian inference with SIR-based compartmental models. Hao et al. [27] proposed an extension of the standard SEIR model, which specifically accounts for asymptomatic and pre-symptomatic individuals. This model, featuring population movement and time varying ascertainment rates, was applied to the initial wave of COVID-19 epidemic in Wuhan, where the confinement measures could be explicitly identified. Bhaduri et al. [28] proposed another modification of the SEIR model, which considers not only asymptomatic infectives, but containment measures and the potential effect of misclassifications (false negative rate) in consequence of testing error. The study demonstrated the effectiveness of the model which was applied to the outbreak in India up to the summer of 2020. As a typical Bayesian approach, the Markov Chain Monte Carlo algorithm was applied for adequate parameter inference in both studies.

With the similar objective for the explanatory analysis of a multi-wave COVID-19 epidemic with change points, some of the related studies are briefly summarized below. Bherwani et al. [29] presented a change point analysis (CPA) that utilizes the Bayesian probability model (BPM) to understand the changing aspects of COVID-19 penetration within the states of India, suggesting that the states which implemented the lockdown before the exponential rise of cases have been able to control the spread of the disease in a very efficient way. Dehning et al. [5] developed a Bayesian framework based on MCMC sampling to the fundamental SIR model for the spread of COVID-19 to characterize central epidemiological parameters, the temporal change of the spreading rate and the potential change points, and the timing and magnitude of intervention effects. Mbuvha and Marwala [30] further applied this approach on a SEIR model, which succeeded in depicting the short-term multi-wave trajectory of South Africa. Other mathematical methodologies have been developed to solve the multi-wave issue as well. Salvadore et al. [31] presented an integro-differential approach to model the wave-like COVID-19 dynamics that accounts for typical evolution parameters. The model started with a Lotka–Volterra equation, and further included a temporal dependency for the transmissibility consisting of the evolving environmental conditions, where vaccinations, virus variants, and various containment strategies were included. The optimal parameters were calculated for each Italian region separately. Beira and Sebastião [32] conceived a comprehensive PSEIRD(S) model to fit the data sets related to the infected, deceased and hospitalized cases reported for the severe post-Christmas outbreak in Portugal. Some of the model parameters fitted vary discretely with time, which indicates different scenarios of various epidemic waves. Cacciapaglia et al. [33] proposed a consistent picture of the wave pattern based on the epidemic renormalization group (eRG) framework. The results demonstrated that the key to control the arrival of the next wave of the pandemic is the strolling period between waves. As a direct extension of [33], Cot et al. [34] employed the eRG framework supplemented with the implementation of vaccination to reproduce and predict the COVID-19 pandemic across the U.S. The study suggested that vaccinations alone are not enough to curb the current and the next waves, and strict social distancing measures should keep being enacted until herd immunity is achieved.

In our study, we aim to solve the issues of explanatory analysis on the long-lasting multi-wave epidemic with the ongoing large-scale vaccination campaign; thus, we integrate Bayesian inference with an extended compartmental SEIR model featuring the new factor of vaccination and the time-varying transmissibility [35] that aligns with multiple waves in the full epidemic trajectory. To interpret the time dependence, we exert specific change points which are potentially driven by assorted causes like governmental containment measures, nonpharmaceutical interventions [36], clustering events, mutation influence or vaccination impact [7]. We propose a Bayesian framework that is mainly designed for the long-term study of the epidemic trajectory featuring multi-waves with several change points. Aiming to yield an explanatory analysis for the dynamic epidemic growth, especially including the recent scenarios with vaccinations, we accomplish a complete process of Bayesian inference of the time-varying transmissibility along with respective uncertainties of epidemiological and statistical parameters based on the real data in Israel. We believe that our modeling framework is simple and applicable to other countries and regions that are suffering from long-term situations of the COVID-19 pandemic.

The structure of this article is organized as follows. Section 2 offers a manifest description of the real data for modeling. Section 3 elaborates the detailed methodologies for our modeling framework. Section 4 demonstrates the complete results for model inference. Conclusions and discussions are given in Section 5.

2. Materials

2.1. Source of Data

The state of Israel is selected as the studied country for our model deployment. We employ the exact number of COVID-19 confirmed positive cases in Israel using the data from 21 February 2020 (the date of the earliest case reported) to 28 September 2021 via the official websites of World Health Organization. The vaccination data are retrieved via Israel Ministry of Health and the dashboard of ’Our World in Data’. All cases are laboratory confirmed.

2.2. Data Description

As the first country in the world to double-vaccinate the majority (over 50%) of its population, Israel is one of the leaders in fighting against the spreading of the pandemic thanks to its prompt roll-out of COVID-19 vaccines [37]. However, the number of newly reported cases still keeps spiraling in the fourth wave of the epidemic, though Israel has tried almost every effort of intervention, including severe extent of social distancing, banning gathering, workplace and school closures, travel controls, contract tracing, and restrictions on internal movement.

The first case of COVID-19 in Israel was officially confirmed on 21 February 2020 when a female citizen tested positive after returning from her quarantine on the highly infected Diamond Princess ship in Japan [38]. The infection remained at a steady pace under mild interventions until a sudden surge occurred in the third week of March. On 19 March, a national state of emergency was officially declared, which brought about legal enforcement of lockdown restrictions [39]: ’staying home unless absolutely essential.’ The lockdown revealed effectiveness when the number of newly infected cases substantially fell into single digits in May, and thus, the government gave the approval of reopening the economy with restrictions gradually being eased [40].

However, the epidemic staged a comeback with the daily number climbing up to 100 again at the end of May, indicating a more influential and persistent second wave of the pandemic. Numerous regional or nationwide lockdowns were carried out aperiodically with restrictions of both governmental interventions and non-pharmaceutical interventions in consideration of the continuous tradeoff between containment and economic reopening, which unexpectedly led to a rapid escalation starting from August. On 10 September 2020, Israel became the country with the highest rate of COVID-19 infections per capita [41]. A new record of over 9000 daily cases was reached at the end of September, and the second wave of epidemic remained unpromising, owing to various reasons, including the summer vacation and several national holidays, large-scale religious events, protests voicing frustration over the response of the local government to the pandemic [42], and the worldwide epidemic explosion (especially in the U.S. and Europe, which are closely connected to Israel by international travel for the recall of Israelis abroad). Confidence decayed as people became impatient with immoderate restrictions and less alert to the spreading of COVID-19 because of its comparatively low severity rate; therefore, the daily cases scarcely declined.

Right after that, a more serious third wave took place from mid-November 2020. On 24 December 2020, the government declared a third nationwide lockdown (that began on 27 December) with more rigorous restrictions, including indoor gatherings of fewer than 10 people personal traveling limited within a 1000-meter radius from home, etc.; see [43]. Aiming to guard against the potential burst from celebrations of Christmas and New Year’s Eve, the lockdown was strictly implemented for 12 days until a set of tightened restrictions was renewed with updates, such as indoor gatherings of 5 people and the complete closure of Ben Gurion International Airport. The tightest lockdown, which consists of short-term curfews at times, was eventually extended to 7 February 2021. However, the interventions hardly effected a radical cure on the resurging of the epidemic, with a new daily record of over 10,000 reported cases. Though costly in economy, the third lockdown was the checkpoint for the majority of Israelis to become vaccinated against COVID-19. After the first dose of Pfizer vaccine administered on 19 December 2020 [8], a tremendous escalation of vaccinations appeared with the support of Green Pass policies implemented by the Israeli government on 21 February 2021, which granted exclusive permission for entry into hotel, theaters, gyms, and religious places, such as synagogues, etc.; see [44]. In order to return to normal life, most Israel residents restored their motivation on being fully vaccinated (i.e., two doses); hence, substantial drops were observed starting from the end of March. From April to June 2021, the period after the nationwide vaccination campaign, there unfolded a stable situation, which manifestly benefited from the high vaccinated rate of the population.

Technically, Israel almost reached herd immunity in June when the daily new cases were down to 0 for some days and people’s lives became almost normal; however, in the last phase of June 2021, daily case numbers began rising again. Mainly caused by the prevalent SARS-CoV-2 B.1.617.2 Delta variant with a much higher epidemiological spreading rate [45], Israel stepped into the fourth wave inevitably [46], and mild interventions with basic restrictions, such as requirement of mask-wearing and travel limitation, were reinstated to slow down the raging transmission of the Delta variant. Up to October 2021, with thousands of daily infected individuals, the ’fourth wave’ has not yet come to an end and the circumstance of future epidemic in Israel remains unclear. Figure 1 briefly demonstrates the multi-wave trajectory of daily new cases in Israel. Even though Israel is the pioneering country on the third vaccination for the booster shot, which practically started in July 2021, the statistics of reinfected individuals and the protection effectiveness remain to be validated. Therefore, in our study, we only focus on the vaccination data for the first doses (over 6 million) and the second doses (over 5.5 million), while excluding the third doses since we temporarily do not consider the complex situation with reinfected individuals.

3. Methods

To technically present the Bayesian framework designed for our epidemic analysis, two statistical modeling parts are elaborated, respectively, in the following context.

3.1. Epidemiological Model

For the analysis of the COVID-19 epidemic, compartmental models are widely used to quantify the transmissibility of disease, as well as depicting the dynamics between infected and uninfected populations. In our study, we apply the susceptible–exposed–infected–recovered (SEIR) compartmental model that has been dominant in the COVID-19 modeling literature. The SEIR model depicts the transition of individuals between 4 stages of conditions, namely, S —susceptible, E—exposed, I—infected, and R—recovered. Note that the difference between compartments of E and I mainly relies on whether the infection of an individual is confirmed. If not, we consider it an exposed individual who remains in incubation and will probably become infectious, while an infected individual should be removed from the compartment of the susceptibles. For the study of COVID-19, the additional compartment of ’E’ is necessary, considering the non-negligible incubation period under high transmissibility of the disease. Hence, we select the SEIR model, rather than the simpler three-compartment susceptible–infectious–recovered (SIR) model, to be the fundamental epidemiological model in our Bayesian framework.

From the perspective of statistics, the SEIR model can be interpreted as a four-state Markov chain model [30]. The dynamic behavior of the SEIR model relies on solving the following system of ordinary differential equations which depict the analytic trajectory of the infectious disease:

$$\begin{array}{cc}\frac{dS}{dt}\hfill & =-\frac{\lambda SI}{N},\hfill \\ \frac{dE}{dt}\hfill & =\frac{\lambda SI}{N}-\sigma E,\hfill \\ \frac{dI}{dt}\hfill & =\sigma E-\mu I,\hfill \\ \frac{dR}{dt}\hfill & =\mu I,\hfill \end{array}$$

(1)

where $S,E,I,$ and R represent the populations of the susceptible, exposed, infected and recovered, respectively. N is the total population, i.e., $N=S+E+I+R$. $\lambda$ and $\mu$ correspond to the transmission rate and the recovery rate, while $\sigma$ is set to be the rate for individuals in incubation becoming infectious. Thus, $\xi =\frac{1}{\sigma }$ and $T_I=\frac{1}{\mu }$ stand for the incubation period [47] and infectious period; see [30]. The basic reproductive number ($R_0$) is defined by $\lambda /\mu$.

Aiming to analyze the multiple waves in the epidemic trajectory of COVID-19, we consider the modeling design by Dehning et al. [5] as a piecewise SEIR model segmented by change points. Different sets of epidemiological parameters are applied for corresponding periods, that is, $\{{\lambda }_i,{\sigma }_i\}$ for the $(i+1)-$th period. However, since no breakthrough has been yielded in the clinical treatment of COVID-19, we simply retain the value of $\mu$ throughout the whole study period.

In order to fit the current development of COVID-19 pandemic, we specially extend the typical SEIR models with new variables of vaccinations. Considering the fact that most types of the COVID-19 vaccine employed worldwide (including Pfizer-BioNTech, Moderna, Astrazeneca, Johnson & Johnson, etc.) require two doses of injection, we add two variables, $V_1$ and $V_2$, to measure the numbers of daily new cases that have completed the first and the second doses of injection, respectively, and ${\alpha }_1$ and ${\alpha }_2$ as the corresponding rates of vaccine efficacy. Note that this extended version of SEIR model differs from the typical compartmental models with vaccination, e.g., SIR-V [9,10,48], since we intend to apply the exact numbers of daily vaccination. Thanks to the open source of precise data on daily vaccinated cases, there is no need to treat V as an additional compartment with unknown values which require to be approximated via the complete dynamic system. Our extended version of SEIR model is summarized as follows:

$$\left\{\begin{array}{cc}\frac{dS}{dt}\hfill & =S\left(t\right)-S(t-1)=-\frac{{\lambda }_i{S}_{t-1}{I}_{t-1}}{N}-{\alpha }_1{V}_{1,t-1}-({\alpha }_2-{\alpha }_1)V_{2,t-1},\hfill \\ \frac{dE}{dt}\hfill & =E\left(t\right)-E(t-1)=\frac{{\lambda }_i{S}_{t-1}{I}_{t-1}}{N}-{\sigma }_i{E}_{t-1},\hfill \\ \frac{dI}{dt}\hfill & =I\left(t\right)-I(t-1)={\sigma }_i{E}_{t-1}-\mu I_{t-1},\hfill \\ \frac{dR}{dt}\hfill & =R\left(t\right)-R(t-1)=\mu I_{t-1}+{\alpha }_1{V}_{1,t-1}+({\alpha }_2-{\alpha }_1)V_{2,t-1},\hfill \end{array}\right.t\in (T_i,T_{i+1}]$$

(2)

where $i=0,1,\cdots ,k$. Correspondingly, the change points are defined as $T_i$, forming $k+1$ periods as $(T_i,T_{i+1}],\mathrm{where}{T}_0=0\mathrm{and}{T}_{k+1}=\mathrm{\Lambda }$ represent the first day and the ending day of our complete study. The full time series of one compartment is an ordered joint set of values in $(k+1)$ periods, e.g., $S=\{S\left(0\right),\cdots ,S\left(T_1\right),S(T_1+1),\cdots ,S\left(T_2\right),\cdots \cdots ,S\left(T_k\right),\cdots ,S\left(\mathrm{\Lambda }\right)\}$. In our study, we solve Equation (2) with a discrete time step $\Delta t=1$ day and $k=4$ for the real data of Israel.

3.2. Model for the Reported Cases

Similar to the definition from a basic SEIR model, the value for the number of daily incidences is proportional to the number of people leaving the E compartment, i.e., $E\left(T\right)-E(T+1)+S\left(T\right)-S(T+1)$; see [14]. Additionally, we set a new parameter $\delta$ for the ratio of reported cases, i.e., $\delta \left(E\left(T\right)-E(T+1)+S\left(T\right)-S(T+1)\right)$. Following the former work of [5], several auxiliary parameters, including the reporting delay D, the amplitude $f_w$ and phase ${\mathrm{\Phi }}_w$ for weekly modulation, are also considered in the modeling process to precisely depict the non-negligible delay and systematic variation of case reporting (i.e., lower numbers at weekends and accumulating reports during the weekdays). Note that the reporting delay is essential since there always exist circumstances where the individuals take rapid tests by themselves without uploading the positive results in time. Hence, the final expression for simulated daily reported cases from our model could be written as follows:

$$C_T={\gamma }_{T-D}(1-f\left(T\right)),\mathrm{for}T>D,$$

(3)

where ${\gamma }_{T-D}=\delta \left(E(T-D)-E(T-D+1)+S(T-D)-S(T-D+1)\right)$ implies the number of incidences for the corresponding day, and $1-f\left(T\right)$ where $f\left(T\right)=\left(1-f_w\right)·\left(1-\left|sin\left(\frac{\pi }{7}T-\frac{1}{2}{\mathrm{\Phi }}_w\right)\right|\right)$ as the reporting fraction. See Figure 2 for a sample plot of $f\left(T\right)$, where a 7-day cyclic period is demonstrated manifestly.

3.3. Bayesian Framework

Mainly following the statistical framework of [5], we construct a typical process of Bayesian inference for the estimation of parameters in our extended SEIR models. The Bayesian framework practically realizes the posterior inference of parameters by updating prior beliefs through customized likelihood [17]. Theoretically, the Bayesian inference is generally based on the Bayes theorem [49] as follows:

$$\mathbb{P}(\theta \mid \mathit{X},M)=\frac{\mathbb{P}(\mathit{X}\mid \theta ,M)\mathbb{P}\left(\theta \right)}{\mathbb{P}\left(\mathit{X}\right)}\propto \mathbb{P}\left(\theta \right)\mathbb{P}(\mathit{X}\mid \theta ,M),$$

(4)

where $\mathbb{P}(\mathit{X}\mid \theta ,M)$ is the posterior of the parameter vector $\theta$, given the model M. $\mathit{X}$ stands for the observed data, which will be denoted as $\widehat{C_T}$ (for the reported cases of day T) in the following context, while $\mathbb{P}\left(\theta \right|\mathit{X},M)$ is the likelihood and $\mathbb{P}\left(\mathit{X}\right)$ the evidence. In our study, $\theta$ consists of 24 parameters,

$$\theta =\{{\lambda }_i,{\sigma }_i,T_k,\mu ,D,I_0,E_0,\eta ,f_w,{\mathrm{\Phi }}_w,{\alpha }_1,{\alpha }_2,\delta \},\mathrm{for}\mathrm{i}=0,1,2,3,4,\mathrm{k}=1,2,3,4.$$

In Bayesian inference, the likelihood is a critical step that partially determines the effectiveness of the modeling. Here, in our Bayesian framework, the likelihood manifests the probability of observing the real data of reported cases. Suggested by Dehning et al. [5], a Student’s T distribution affiliated to the daily new cases is used to construct the likelihood function:

$$\begin{array}{c}\mathbb{P}\left({\widehat{C}}_T\mid \theta \right)\sim \mathrm{Student}-\mathrm{T}\left(\nu =4,\mu =C_T\left(\theta \right),\sigma =\eta \sqrt{C_T\left(\theta \right)}\right),\end{array}$$

(5)

where $\mu$ is the location parameter, $\sigma$ is the scale parameter, and $\nu$ is the degree of freedom. Note that the parameter $\eta$, which is also in the parameter space $\theta$, is a scaling factor tuning the spread of the distribution. The Bayesian inference based on the MCMC algorithm provides parameter estimation that lessens the discrepancy between the fitted number $C_T\left(\theta \right)$ (for any time point T) and the corresponding real data ${\widehat{C}}_T$; thus, it statistically quantifies the similarity between the predicted cases derived from the model (3) and the real-world time series for observed reported cases. The distribution around the mean is similar to a Gaussian one; however, with heavy tails, it enhances the robustness of MCMC with respect to outliers, i.e., reporting noise [50]. The spread (i.e., the scale parameter $\sigma$), motivated by observation noise through subsampling [51], is highly dependent on the number of daily cases. This design forms a variance proportional to the mean of the Student’s T distribution. Note that the general assumption of neglecting noises in the dynamic process is used since the SEIR compartmental model is fully deterministic.

3.4. Prior Specification

The empirical prior knowledge is taken into consideration in a Bayesian setting while obtaining reasonable parameter estimation for epidemic analysis and short-term forecasting. The elicitation of prior knowledge for the epidemiological parameters (e.g., $R_0$), the transmissibility and recovery rate of the disease, the mechanisms of the spreading trajectory, etc., is based on former studies (or expertise knowledge) with preliminary estimation from other countries and regions or from another closely related pandemic contagion. Based on expertise and empirical knowledge, the setting of prior distributions of $\theta$ plays an important role for parameter estimation. On the other hand, for the epidemiological parameters in our model, such as spreading rate $\lambda$, recovery rate $\mu$, infectious rate $\sigma$, change points based on epidemiological periods, etc., their priors are specified with expectation values that are referred from related studies.

The detailed prior distributions of all the parameters are fully listed in

Table 1. Priors on the model parameters.

Parameter	Interpretation	Prior
${\lambda }_0$	Initial transmission rate	LogNormal(log(0.4), 0.5)
${\lambda }_1$	Transmission rate after the 1st change point	LogNormal(log(0.4), 0.5)
${\lambda }_2$	Transmission rate after the 2nd change point	LogNormal(log(0.4), 0.5)
${\lambda }_3$	Transmission rate after the 3rd change point	LogNormal(log(0.4), 0.5)
${\lambda }_4$	Transmission rate after the 4th change point	LogNormal(log(0.4), 0.5)
${\sigma }_0$	Initial infectious rate	LogNormal(log(1/5), 0.5)
${\sigma }_1$	Infectious rate after the 1st change point	LogNormal(log(1/5), 0.5)
${\sigma }_2$	Infectious rate after the 2nd change point	LogNormal(log(1/5), 0.5)
${\sigma }_3$	Infectious rate after the 3rd change point	LogNormal(log(1/5), 0.5)
${\sigma }_4$	Infectious rate after the 4th change point	LogNormal(log(1/5), 0.5)
$\mu$	Recovery rate	LogNormal(log(1/8), 0.2)
D	Reporting delay	LogNormal(log(8), 0.2)
$I_0$	Initial number of infected	HalfCauchy(100)
$E_0$	Initial number of exposed	HalfCauchy(100)
$T_1$	The 1st change point	Normal(7 May 2020, 3)
$T_2$	The 2nd change point	Normal(17 Aug 2020, 3)
$T_3$	The 3rd change point	Normal(31 Oct 2020, 3)
$T_4$	The 4th change point	Normal(20 Jun 2021, 3)
$f_w$	Weekly modulation amplitude	Beta(4.39, 1.88)
${\mathrm{\Phi }}_w$	Weekly modulation phase	VonMises(0, 0.01)
$\delta$	The proportion of cases reported	Beta(1, 2)
$\eta$	Scale factor	HalfCauchy(10)
${\alpha }_1$	Effective rate of the first dose	Normal(0.52,0.01)
${\alpha }_2$	Effective rate of the second dose	Normal(0.95,0.01)

. Mainly following the former study by [5], we select log-normal distributions for the transmission rate $\lambda$ and recovery rate $\mu$ such that the initial basic reproduction number $R_0$ [20,52] equals a median of $3.2$. The infectious rate $\sigma$ is set to have the log-normal prior, which is equivalent to the mean incubation period of 5 days. Additionally, we also choose a log-normal prior for the reporting delay $\mathrm{D}\sim \mathrm{LogNormal}(log(8),0.2)$ to incorporate both the mean value of the incubation period with median 5 days plus an assumed 3 days delay from being infected to being reported; see [53].

We recognize the change points in consideration of the Stringency Index proposed by the research project of ’Oxford COVID-19 Government Response Tracker’ along with checkpoints for the change in governmental intervention policies and nationwide events; see Figure 3. The Stringency Index is a composite measure based on 9 response indicators, including school closures, workplace closures, and travel bans, rescaled to a score value from 0 to 100, which is proportional to the strictness of governmental interventions; see [54].

As a result, we roughly treat 7 May 2020 as the first change point when the exit strategy was implemented after the first wave of epidemic boost; see [40]. From Figure 1, we notice the obvious two-stage pattern of escalation within the duration of the second wave, manifesting that an additional change point should not to be ignored. Indicated by the following trough in the curve of the Stringency Index, 17 August 2020 is picked as the second change point. In the same manner, by recognizing the substantial bottom of a certain period in the trajectory of the Stringency Index, we set 31 October 2020 as the starting point of the durable third wave and 15 June 2021 as the beginning of the current wave (the ’fourth wave’) as well as the division between the pre-Delta and post-Delta periods. Normal distributions with a standard deviation of 3 days are applied to construct the priors for $T_1$, $T_2$, $T_3$ and $T_4$. Note that we keep the same prior for ${\lambda }_0$, ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ and ${\lambda }_4$ (also similar to ${\sigma }_i$) since we are uninformed of the exact transmission rates after the initial phase of the pandemic, which will lead to a greater weight on the data-driven likelihood [30].

The parameters of initial infected individuals $I_0$, initial exposed individuals $E_0$ and the scale factor for the likelihood $\eta$ are constrained by non-informative priors. The half-Cauchy distributions are employed with corresponding expectations.

Similar to the modeling design in [5], we construct a non-negative sine function (applying absolute values) with a total period of 7 days to modulate the number of daily cases; hence, we select a weakly informative prior $\mathrm{Beta}(4.39,1.88)$ (equivalent to the mean of 0.7 and the standard deviation of 0.17) for the amplitude of weekly modulation $f_w$, along with a prior as a nearly flat circular distribution $\mathrm{VonMises}(0,0.01)$ for the phase ${\mathrm{\Phi }}_w$.

We give the proportion of reported cases $\delta$ a weakly-informative $\mathrm{Beta}(1,2)$ prior, which indicates that we are quite uncertain about this parameter, except that it stays within the interval of $(0,1)$ but should not be close to 1. For the remaining two parameters, ${\alpha }_1$ and ${\alpha }_2$, which respectively stand for the efficacy of the first and the second dose of COVID-19 vaccine, we propose normal distributions with corresponding mean values for these priors. According to Pfizer data published in December 2020 [11], the Pfizer-BioNTech vaccine is roughly 52% effective after the first dose with 95% protection rate for fully vaccinated people. Though the Pfizer-BioNTech vaccine is not the only one distributed in Israel (Moderna and Pfizer-BioNTech are the two manufacturers available to the population of Israel), we still adopt these authenticated data since the two vaccines were similarly developed based on the same mRNA sequence, resulting in little difference on clinical efficacy. Table 1 summarizes the prior distributions of all the parameters described above.

3.5. Markov Chain Monte Carlo (Mcmc) with No-U-Turn Sampler (NUTS)

We adopt the MCMC algorithm as a straightforward and plausible approach for the effective inference of the posterior distributions on the parameters constructed; see [55,56]. MCMC works by drawing a sequence of dependent samples from the posterior distribution under circumstances, where it is unlikely for a direct draw of independent samples. With the algorithm, a Markov chain $\{{\theta }_t,t=1,2,\cdots \}$ is created (representing a draw at the state of fictitious time t) such that

$$\mathbb{P}\left({\theta }_{t+1}\mid {\theta }_1,\dots ,{\theta }_t\right)=\mathbb{P}\left({\theta }_{t+1}\mid {\theta }_t\right).$$

(6)

The main advantage of the MCMC approach for our complete model is that the full covariance matrix of $\theta$ is used; see [57]. Posterior inference on $\theta$ is carried out by random samples from the conditional distribution given the observation data, i.e., the reported cases ${\widehat{C}}_T$ of day T, and it follows from Equation (4) that

$$\mathbb{P}(\theta \mid {\widehat{C}}_T)\propto \mathbb{P}({\widehat{C}}_T\mid \theta )\mathbb{P}\left(\theta \right).$$

(7)

As one of the most frequently used methods in practice, the Metropolis–Hastings (MH) algorithm [58,59] is a simple and well-known MCMC method for generating a Markov Chain which converges to the target distribution. The MH generates a candidate ${\theta }_{t\ast }$ from a proposal distribution $\mathbb{P}\left({\theta }_{t\ast }\mid {\theta }_{t-1}\right)$ with acceptance ratio

$$\mathbb{P}\left(accept\left({\theta }_{t\ast }\right)\right)=min\left(1,\frac{\mathbb{P}\left({\theta }_{t\ast }\mid {\widehat{C}}_T\right)\mathbb{P}\left({\theta }_{t\ast }\mid {\theta }_{t-1}\right)}{\mathbb{P}\left({\theta }_{t-1}\mid {\widehat{C}}_T\right)\mathbb{P}\left({\theta }_{t-1}\mid {\theta }_{t\ast }\right)}\right)$$

(8)

In most cases, a symmetric random walk is applied as the common proposal distribution, which is obtained by adding a Gaussian noise to a former state of the parameter. However, the proposal distribution with the random walk behavior often brings about a relatively low sample acceptance and limited sizes of effective samples. To overcome the shortage, the Hamiltonian Monte Carlo algorithm (HMC) was proposed by [60] to reduce random walk behavior by adding an auxiliary momentum variable $r_i$ to each parameter space ${\theta }_i$. The basic idea of HMC is to employ the current gradient to assign a direction heading for the next state with a high probability. In general form, the HMC dynamic system, consisting of model parameters $\theta$ and the auxiliary momentum variables r, is denoted by the Hamiltonian $H(\theta ,r)$, written as follows:

$$H(\theta ,r)=U\left(\theta \right)+K\left(r\right),$$

(9)

where $U\left(\theta \right)=logp(\theta ,r)$ signifies the negative log-likelihood of the posterior. From the perspective of the Hamiltonian, this equation could be referred to as a combination of the potential energy $U\left(\theta \right)$ and kinetic energy $K\left(r\right)$, which is specified by a Gaussian kernel with the covariance matrix $\Sigma$ [61].

$$K\left(r\right)=\frac{r^T{\Sigma }^{-1}r}{2},$$

(10)

and the Hamiltonian dynamics for Equation (10) (at a fictitious time t) is defined as follows:

$$\begin{array}{c}\frac{\partial H}{\partial r_i}=\frac{\partial {\theta }_i}{\partial t},\\ \frac{\partial H}{\partial {\theta }_i}=-\frac{\partial r_i}{\partial t}.\end{array}$$

(11)

During the standard procedure of HMC, the leapfrog integrator discretizes the dynamical trajectory of Equation (10). Within the leapfrog integrator, a doubling process, which includes a half step in the direction of r, followed by a full step (with a discretization stepsize $ϵ$) in the direction of $\theta$, then ending with another half step in the same direction as the first half step, is operated to reach the next point of the path. Following the process of the Metropolis–Hastings algorithm, the newly derived parameters ${\theta }^{\ast }$ take the probability of acceptance as follows [62]:

$$\mathbb{P}\left(\mathrm{accept}\left({\theta }^{\ast }\right)\right)=min\left(1,\frac{\mathbb{P}\left({\theta }^{\ast }\mid {\widehat{C}}_T,H(\theta ,r)\right)}{\mathbb{P}\left(\theta \mid {\widehat{C}}_T,H(\theta ,r)\right)}\right).$$

(12)

The effective performance of HMC strongly relies on selecting adequate values for the discretization stepsize $ϵ$ and the trajectory length L. Reasonable tuning for these two parameters is crucial for efficient sampling. Therefore, the No-U-Turn Sampler (NUTS) is applied as an extension of HMC that eliminates the need for specifying a value of L as well as retaining the capacity to inhibit random walk behavior; see [30,63]. In short, NUTS automates the tuning of $ϵ$ and L. Specifically, $ϵ$ is tuned via an initial burn-in phase by assigning a target mean acceptance probability, while L is chosen by adding extra steps iteratively until the trajectory starts to retrace itself, or descriptively speaking, to make a ’U-turn’.

NUTS suggests a recursive algorithm, which is a revised version of the doubling procedure based on slice sampling [64], thus a slice variable u is introduced with a conditional uniform distribution:

$$P(u\mid \theta ,r)=\left\{\begin{array}{cc}\frac{1}{Z}\hfill & 0\le u\le exp\left\{\mathcal{L}\left(\theta \right)-\frac{1}{2}r·r\right\},\hfill \\ 0\hfill & \mathrm{Otherwise},\hfill \end{array}\right.$$

(13)

where the normalization constant $Z=\int exp\left\{\mathcal{L}\left(\theta \right)-\frac{1}{2}r·r\right\}d\theta$, with $\mathcal{L}$ as the logarithm of the density of $\theta$. This procedure aims to sample above the slice of the posterior within a targeted bandwidth so as to guarantee the acceptance. After resampling, the Leapfrog integrator is implemented with steps taken back and forth to portray the path, which includes the first 1 step, then 2 steps, and 4 steps afterwards, etc. Hence, a balanced binary tree is constructed with leaf nodes corresponding to the position-momentum states. The signal for a stop is when the trajectory from the leftmost to the rightmost nodes of any subtree begins to double back (i.e., make a ’U-turn’). NUTS keeps doubling the tree at a height j until the condition is achieved that all the $2^j-1$ balanced binary subtrees have positive values of height, and the doubling process is halted when one of these subtrees with states associated with the leftmost and rightmost leaves ${\theta }_{-},r_{-}$ and ${\theta }_{+},r_{+}$ satisfies the following condition:

$$\left({\theta }^{+}-{\theta }^{-}\right)·r^{-}<0\mathrm{or}\left({\theta }^{+}-{\theta }^{-}\right)·r^{+}<0.$$

(14)

An upper boundary ${\Delta }_{max}$ (usually, we take a large positive number, such as 1000) is set to be the second rule: we simply stop doubling if the tree includes a leaf node whose state $(\theta ,r)$ satisfies

$$\mathcal{L}\left(\theta \right)-\frac{1}{2}r·r-logu<-{\Delta }_{max},$$

(15)

To summarize, Algorithm 1 elaborates the complete process of our Bayesian modeling framework. Note that for simplicity, Algorithm 1 is written in a pseudocode in Visual Basic.

Algorithm 1 Complete modeling framework

i. Sample from parameter priors set in Table 1 to obtain $\theta$. Given initialization of $S_0=N-I_0-E_0,R_{int}=0$ with values of N and $\mathrm{\Lambda }=T_{k+1}$ derived from the applied real data. for$t\in [T_i+1,T_{i+1}]$, $i\leftarrow 0$ to k do: $\left\{\begin{array}{cc}\hfill \frac{dS}{dt}& =S\left(t\right)-S(t-1)=-\frac{{\lambda }_i{S}_{t-1}{I}_{t-1}}{N}-{\alpha }_1{V}_{1,t-1}-({\alpha }_2-{\alpha }_1)V_{2,t-1},\hfill \\ \hfill \frac{dE}{dt}& =E\left(t\right)-E(t-1)=\frac{{\lambda }_i{S}_{t-1}{I}_{t-1}}{N}-{\sigma }_i{E}_{t-1},\hfill \\ \hfill \frac{dI}{dt}& =I\left(t\right)-I(t-1)={\sigma }_i{E}_{t-1}-\mu I_{t-1},\hfill \\ \hfill \frac{dR}{dt}& =R\left(t\right)-R(t-1)=\mu I_{t-1}+{\alpha }_1{V}_{1,t-1}+({\alpha }_2-{\alpha }_1)V_{2,t-1},\hfill \end{array}\right.$ end for return$S=\left\{S\right(1),\cdots ,S(\mathrm{\Lambda }\left)\right\}$, $E=\left\{E\right(1),\cdots ,E(\mathrm{\Lambda }\left)\right\}$, $I=\left\{I\right(1),\cdots ,I(\mathrm{\Lambda }\left)\right\}$, $R=\left\{R\right(1),\cdots ,R(\mathrm{\Lambda }\left)\right\}$ ii. Calculate reported cases $C_T\left(\theta \right)$ via Equation (3):       $C_T\left(\theta \right)=\left\{\begin{array}{ccc}& {\gamma }_T,\hfill & \hfill \mathrm{for}T\in \{1,2,\cdots ,D\}\\ & {\gamma }_{T-D}(1-f\left(T\right)),\hfill & \hfill \mathrm{for}T>D\end{array}\right.$ iii. Perform MCMC-NUTS with the Student-T likelihood of Equation (5). Set $\mathcal{L}\left(\theta \right)=log\mathbb{P}(\theta \mid {\widehat{C}}_T)$ where $\mathbb{P}(\theta \mid {\widehat{C}}_T)\propto \mathbb{P}({\widehat{C}}_T\mid \theta )\mathbb{P}(\theta )$. Given non-informative initialization of $ϵ,L$ and conduct the following iteration steps: for$m\leftarrow 1$ to L do:       1. $r_0\sim \mathcal{N}(0,I)$       2. $u\sim \mathrm{Uniform}\left([0,exp\{\mathcal{L}\left({\theta }_{m-1}\right)-\frac{1}{2}{r}_0·r_0)\}\right])$       3. $r_{-}=r_0=r_{+},n=1,s=1,{\theta }_{-}\leftarrow {\theta }_{m-1},{\theta }_{+}\leftarrow {\theta }_{m-1},{\theta }_m\leftarrow {\theta }_{m-1}$    while s = 1:       Choose a direction $v_j\sim \mathrm{Uniform}\left(\{-1,1\}\right\})$       if $v_j=-1$:          ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime },n^{\prime },s^{\prime }\leftarrow$ BuildTree$({\theta }_{-},r_{-},u,v_j,j,ϵ)$       else:          ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime },n^{\prime },s^{\prime }\leftarrow$ BuildTree$({\theta }_{+},r_{+},u,v_j,j,ϵ)$       end if       if $s^{\prime }=1$ then          Set ${\theta }_m\leftarrow {\theta }^{\prime }$ with acceptance probability $min\{1,\frac{n^{\prime }}{n}\}$       end if       1. $n\leftarrow n+n^{\prime }$       2. $s\leftarrow s^{\prime }\mathbb{I}\left[\left({\theta }^{+}-{\theta }^{-}\right)·r^{-}\ge 0\right]\mathbb{I}\left[\left({\theta }^{+}-{\theta }^{-}\right)·r^{+}\ge 0\right]$       3. $j\leftarrow j+1$    end while end for return$\widehat{\theta }={\theta }_L$. function BuildTree$(\theta ,r,u,v,j,ϵ)$ if$j=0$then    One Leapfrog step in the direction v.    1. ${\theta }^{\prime },r^{\prime }\leftarrow Leapfrog(\theta ,r,vϵ)$    2. $n^{\prime }\leftarrow \mathbb{I}\left[u\le exp\left\{\mathcal{L}\left({\theta }^{\prime }\right)-\frac{1}{2}{r}^{\prime }·r^{\prime }\right\}\right]$    3. $s^{\prime }\leftarrow \mathbb{I}\left[u<exp\left\{{\Delta }_{max}+\mathcal{L}\left({\theta }^{\prime }\right)-\frac{1}{2}{r}^{\prime }·r^{\prime }\right\}\right]$    return ${\theta }^{\prime },r^{\prime },n^{\prime },s^{\prime }$. else    Build the left and right subtrees: ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime },n^{\prime },s^{\prime }\leftarrow$ BuildTree$(\theta ,r,u,v,j-1,ϵ)$    if $s^{\prime }=1$ then       if v = −1 then          ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime \prime },n^{\prime \prime },s^{\prime \prime }\leftarrow$ BuildTree$({\theta }_{-},r_{-},u,v,j-1,ϵ)$       else          ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime \prime },n^{\prime \prime },s^{\prime \prime }\leftarrow$ BuildTree$({\theta }_{+},r_{+},u,v,j-1,ϵ)$       end if       1. Set ${\theta }^{\prime }\leftarrow {\theta }^{\prime \prime }$ with acceptance probability $\frac{n^{\prime \prime }}{n^{\prime }+n^{\prime \prime }}$       2. $s^{\prime }\leftarrow s^{\prime \prime }\mathbb{I}\left[\left({\theta }^{+}-{\theta }^{-}\right)·r^{-}\ge 0\right]\mathbb{I}\left[\left({\theta }^{+}-{\theta }^{-}\right)·r^{+}\ge 0\right]$       3. $n^{\prime }\leftarrow n^{\prime }+n^{\prime \prime }$    end if    return ${\theta }_{-},{\theta }_{+},r_{-},r_{+},{\theta }^{\prime },n^{\prime },s^{\prime }$ end if function Leapfrog$(\theta ,r,ϵ)$       1. $\tilde{r}\leftarrow r+(ϵ/2){\nabla }_{\theta }\mathcal{L}\left(\theta \right)$       2. $\tilde{\theta }\leftarrow \theta +ϵ\tilde{r}$       3. $\tilde{r}\leftarrow \tilde{r}+(ϵ/2){\nabla }_{\theta }\mathcal{L}\left(\tilde{\theta }\right)$ return$\tilde{\theta },\tilde{r}$ iv. Derive $C_T(\widehat{\theta })$ for fitted numbers and $C_{T^{\prime }}(\widehat{\theta })$ ($T^{\prime }$ as a short duration right after the study period T) for the prediction.

4. Results

We propose a novel Bayesian framework combining a modified piecewise SEIR model with an empirical model on daily reported and vaccinated cases in Israel, which captures both the change points in the long-term multi-wave trajectory and oscillations in the short-term period. Reporting delay and reporting fraction are considered to justify the fitness of the model.

To approximate the posterior distributions, four independent Markov Chains are performed with 1000 burn-in (tuning) steps and 5000 samples for each chain of the sampler. In our study, we apply the model and algorithm to the real data from 21 February 2020 to 28 August 2021 (as the ending day $\mathrm{\Lambda }$), while leaving 30 days (from 29 August 2021 to 27 September 2021) for prediction. According to the population of Israel in 2020, we use an approximate number of 9.2 million as the known value for N.

The detailed posterior results of our extended SEIR model using NUTS are presented as follows. Figure 4 portrays the prior and posterior distributions of our modified SEIR model parameters. Figure 5 demonstrates the comparison between the real multi-wave trajectory and the results of daily incidence ${\gamma }_T=\delta (E\left(T\right)-E(T+1)+S\left(T\right)-S(T+1))$ obtained directly from our modified SEIR model, where a smooth curve (free from noise) is fitted without the calibration of the reporting delay and the weekly modulation (which definitely adds noise to the yielded SEIR results).

Table 2. Fitting results for model parameters and derived epidemiological parameters.

Parameter	Mean	Standard Deviation	95% Confidence Interval
${\lambda }_0$	0.62	0.07	[0.50, 0.76]
${\lambda }_1$	0.34	0.04	[0.28, 0.45]
${\lambda }_2$	0.36	0.05	[0.25, 0.45]
${\lambda }_3$	0.49	0.07	[0.37, 0.64]
${\lambda }_4$	1.02	0.13	[0.82, 1.31]
${\sigma }_0$	1.05	0.42	[0.47, 2.09]
${\sigma }_1$	0.32	0.11	[0.16, 0.59]
${\sigma }_2$	1.83	0.63	[0.94, 3.35]
${\sigma }_3$	0.33	0.04	[0.25, 0.42]
${\sigma }_4$	1.26	0.49	[0.57, 2.45]
$\mu$	0.31	0.05	[0.23, 0.41]
D	7.68	0.92	[5.72, 9.30]
$I_0$	1.51	1.87	[0.05, 6.45]
$E_0$	1.59	2.02	[0.04, 6.94]
$T_1$	8 May 2020	0.94	[6 May 2020, 10 May 2020]
$T_2$	20 August 2020	0.68	[19 August 2020, 22 August 2020]
$T_3$	1 November 2020	1.51	[28 October 2020, 4 November 2020]
$T_4$	16 June 2021	1.63	[13 June 2021, 20 June 2021]
$f_w$	0.58	0.06	[0.47, 0.69]
${\mathrm{\Phi }}_w$	3.71	2.72	[1.74, 8.60]
$\delta$	0.20	0.01	[0.18, 0.22]
$\eta$	19.79	0.89	[18.11, 21.64]
${\alpha }_1$	0.52	0.03	[0.47, 0.57]
${\alpha }_2$	0.95	0.01	[0.93, 0.97]
${\left(R_0\right)}_1$	2.03	0.24	[1.64, 2.56]
${\left(R_0\right)}_2$	1.10	0.01	[1.07, 1.13]
${\left(R_0\right)}_3$	1.17	0.04	[1.11, 1.24]
${\left(R_0\right)}_4$	1.59	0.04	[1.50, 1.68]
${\left(R_0\right)}_5$	3.33	0.32	[2.79, 4.03]
${\xi }_1$	1.10	0.43	[0.48, 2.12]
${\xi }_2$	3.50	1.14	[1.70, 6.12]
${\xi }_3$	0.61	0.20	[0.30, 1.07]
${\xi }_4$	3.11	0.40	[2.38, 3.93]
${\xi }_5$	0.91	0.35	[0.41, 1.75]

shows the summary of the final fitting results, including the transmission rates ${\lambda }_i,i=0,1,2,3,4$ and the infectious rates ${\sigma }_i,i=0,1,2,3,4$ in five periods, along with the overall recovery rate $\mu$, the proportion of cases reported $\delta$, and the reporting delay D. The corresponding basic reproduction numbers ${\left(R_0\right)}_i,i=1,2,3,4,5$ and the incubation periods ${\xi }_i,i=1,2,3,4,5$ are also derived from the obtained parameters. From the results, we notice that the transmission rates vary significantly from ${\lambda }_3$ to ${\lambda }_4$, as well as the similar trend that the infectious rates surge drastically from ${\sigma }_3$ to ${\sigma }_4$. The fact reflects that $T_4$ is a remarkable change point for the transmissibility of COVID-19, which denotes the actual inflection point of the epidemic period before and after the diffusion of the Delta variant. In a same manner, we could further infer that there exists a substantial difference for the values of $R_0$ and $\xi$ of the pre-Delta and post-Delta periods. With higher transmissibility, a larger reproduction number and the lower incubation period are estimated for the spreading of the Delta variant in Israel.

Figure 6 shows the projections of daily new cases throughout the complete studied period based on our extended SEIR models with the four change points. Figure 7 shows the corresponding curves for cumulative cases. Meanwhile, the mean values of change points are inferred as $T_1=$ 8 May 2020, $T_2=$ 20 August 2020, $T_3=$ 1 November 2020 and $T_4=$ 16 June 2021. These inferred dates match the actual timing of the restrictions issued by the government at different stages, including reopening malls and markets, resuming outdoor cultural performances and religious gatherings, school returning and even the drop of the indoor mask mandate; see [40,65,66,67], respectively. Synchronicity partially reveals the connection between change points and the timing of easing nonpharmaceutical interventions, which is consequently followed by new waves of epidemic outbreak.

The empirical vaccination data mainly promote the goodness of fit in the fourth wave, in which the majority of Israel residents have been injected with at least one dose of vaccine, resulting in an increasing scale of immunity. Validated by an adequate goodness of fit (with R-squared = 0.71 and the p-value of ${\chi }^2$-test as $0.26>0.05$), our 30-day prediction demonstrates its success in depicting the fourth wave and the tendency of rapid subsidence.

In order to reveal the certain novelty added by the vaccination effect (i.e., $V_1\mathrm{and}{V}_2$ along with coefficients ${\alpha }_1\mathrm{and}{\alpha }_2$) featured in our modified SEIR model, we additionally perform a numerical experiment that simulates on a simpler model (with the same initialization of compartments) with no application of the empirical vaccination data. Without loss of generality, we only examine the trajectories in the latest period (i.e., $[T_4,\mathrm{\Lambda }]$), where a substantial proportion of Israelis has been vaccinated for at least one dose, which is expected to show the most remarkable difference throughout the full study period. Figure 8 portrays the comparisons between two models, and the fitted results clarify that the simpler model (without vaccination) might give rise to an overestimation of daily new cases as well as the transmissibility, while the epidemic peak is predicted over one week later than the actual timing. This further validates the advantage of our modified SEIR model featuring the vaccination effect.

5. Discussion

We utilized a Bayesian inference framework for modeling time-varying transmission rates of COVID-19 in Israel. As is shown in the results, the time-varying transmission rates offer explanatory analysis of multiple waves of the epidemic growth, along with assessing the impact of governmental nonpharmaceutical interventions and the effectiveness of vaccination in practice. It follows that the spreading of COVID-19 in the first period (from 21 February 2020 to $T_1$) has comparatively high transmissibility, though the epidemic growth in Israel was under effective control of the nationwide lockdown. However, owing to capricious intervention policies balancing between containment and economic reopening, in the second and the third periods (from $T_1$ to $T_2$, and from $T_2$ to $T_3$, respectively), the country encountered two-stages of epidemic surge with similar mild transmission rates but widely varied incubation rates, which led to different acceleration of daily new cases. The third wave of epidemic in Israel (equivalent to the fourth studied period from $T_3$ to $T_4$) was observed to have a relatively higher transmission rate and a long incubation period, resulting in the long-lasting growth curve containing a complete process of epidemic surge and subsidence influenced by the continuous support of flexible mitigation measures. Nevertheless, principally in consequence of the raging Delta variant, the fourth wave (equivalent to the last studied period from $T_4$ to the ending day) kept testifying a tremendous soaring with the highest transmissibility. Though the massive vaccination campaign was expected to help Israelis ’return to normal’ shortly while the exit strategies were gradually performed, the governmental interventions were almost lifted until the comeback pandemic, which is dominant with the Delta variant. However, thanks to the vaccination process, our prediction forecasts a rapid decrease right after the crest around early September 2021, which corresponds to the actual circumstances.

Compared to the related studies with similar aim, our model presents both simplicity and functionality. Among them, the PSEIRD(S) model given by Beira and Sebastião [32] revealed its capability to explain the capricious epidemic wave. The model fitting was based on the data sets of daily infected, daily deceased, hospitalized and ICU patients in Portugal. However, considering the complexity of this modified compartmental model, the parameter estimation, which is implemented by the exhaustive Runge–Kutta method, requires high level of computational consumption. Though the process of vaccination was simulated in the numerical study, the real data of daily vaccinated cases has not been applied in the model. The epidemic renormalization group (eRG) framework employed in the work of Cacciapaglia et al. [33] turned out to be a useful tool for the prediction of the next wave. As a modification performed by Cot et al. [34], the vaccination was also utilized by a random fraction $R_v$ of nodes set to the state V. This approach is based on the analysis of the time evolution of the number of total cases and the symmetries revealed by the epidemic trajectory in order to fetch information from the data on the specific conditions of each studied country. Unlike compartmental models, this mathematical framework does not account for the exponential growth and could not be applied directly to the first two waves due to the uncertainty in their relative normalization. It is simple in form, but not as straightforward as compartmental models in terms of explainability, and the computational consumption required for the numerical optimization that deals with a great amount of parameters is also non-negligible. Salvadore et al. [31] presented a novel model, which is based on the Lotka–Volterra equation and takes into account the environmental factors, including vaccinations, virus variants, and different containment measures, as an integro-differential approach to portray the wave-like COVID-19 dynamics. The optimization algorithm for parameter fitting is differential evolution, a metaheuristic method performed by an exhaustive search for the global minima. Despite the functionality, this model still lacks computational efficiency to some extent, and takes no practical usage of the empirical vaccination data. For our method, we mainly modified the Bayesian framework on the SEIR compartmental model by piecewise modeling (divided by change points) and the additional using of vaccination data. The basic framework was originally proposed by Dehning et al. [5] and was further applied by Mbuvha and Marwala [30] with the extension on SEIR model; however, their previous studies only validated the effectiveness on short-term trajectories of Germany and South Africa before the vaccination campaign had taken place. The Bayesian framework provides efficient calculation of parameter inference and the practical usage of prior information, which takes obvious advantage over other methods in the related studies mentioned above. Moreover, there is no explicit research on the explanatory analysis of the complete multi-wave epidemic in Israel. Hence, the novelty of our study appears as an application on the full-time COVID-19 trajectory of a typical country with multi-wave pattern and high vaccinated rate of the population.

Although our model has derived an explanatory analysis of the multi-wave epidemic trajectory, some technical points remain to be further discussed. Firstly, despite the good fitness for the long-term history data (from 21 February 2020 to 28 August 2021), our model shows an relative underestimation of the short-term prediction (from 29 August 2021 to 27 September 2021). A direct explanation is that we overestimated the deduction of susceptible individuals caused by the increasing population of vaccinated residents since there exists an unknown dynamic of second infection beyond the validity of vaccine protection and waning of immunity. Based on the current SEIR compartmental model, a modification could be proposed with an additional dynamic between compartments of ‘R’ and ‘S’. As a reference for studying the epidemiological dynamics with reinfection, Malkov [68] proposed a susceptible–exposed–infectious–resistant–susceptible (SEIRS) model with three different ways of modeling reinfection, including the model with a constant immunity waning rate, the infected one with a milder form of the disease, and the individuals that are resistant at some date who become susceptible again. Simulations were presented to explore the difference in the spreading of COVID-19 under the reinfection and no-reinfection scenarios as well as the effectiveness of the mitigation measures. The study argued that the dynamics of the reinfection and no-reinfection scenarios are indistinguishable before the epidemic prime. The models also prevailed in that the mitigation measures delay not only the peak of infection, but the moment when a remarkable difference between the reinfection and no-reinfection scenarios emerges as well. Results showed the robustness to various modeling assumptions (e.g., different waning rates, and various basic reproduction numbers).

Secondly, although we showed the smoothness of the results from our modified SEIR model (Figure 5), oscillations are universally observed throughout the final fitted numbers of daily new cases; see Figure 6. The reason behind this is that the assumption of $C_T$ (Equation (3)) is constructed by ${\gamma }_T$ with the corresponding calibration of the reporting delay D and the reporting fraction $f\left(T\right)$ that consists of the weekly modulation $f_w,{\mathrm{\Phi }}_w$. This setting diminishes the smoothness of the results of ${\gamma }_T$ that is directly yielded from our modified SEIR model, but it facilitates the effective fitting to the noised real-world data via Bayesian inference.

Thirdly, in the collection of real data, there exist a few outliers which turn out to be the noise in the modeling process, e.g., 0 case for the days (6 September 2021 and 7 September 2021) when there were no official statistics recorded, and an extremely high number of cases for the following day (8 September 2021) as the calibration. In our study, we have no preprocessing for the limited outliers. However, the robustness has not been verified in the circumstances with a larger amount of outliers, which is commonplace in statistics.

In general, the limitation of our modeling framework is mainly about the range of functions. It is true that we are seeking the balance between simplicity and utility; thus, we only considered a typical four-compartmental SEIR model for depicting the multi-wave trajectory. Nevertheless, particularly for the pandemic of COVID-19, it is somehow worthwhile to extend the current dynamic system with additional compartments, such as asymptomatic individuals and hospitalized patients. For a high-transmissible contagion with a long incubation period, the populations of asymptomatic and hospitalized individuals, along with time-varying containment strategies and the reinfection scenarios, are important characteristics required to be studied by a more complex model, but we believe that these points of interest could be easily realized by extending the dynamic system in our modeling framework. Moreover, indicated by a recent report [69], there exists some difference in vaccine efficacy for new variants (the SARS-CoV-2 B.1.617.2 Delta variant); thus, another limitation of our model is that it has not taken into account for the changing vaccine. Following the idea of piecewise modeling, a straightforward modification is to assign the coefficients of vaccine efficacy (namely ${\alpha }_1$ and ${\alpha }_2$) with distinct values for study periods before and after the corresponding change point. However, this remains a simple but rough solution, which only considers fixed values for pre-Delta and post-Delta vaccine efficacy, while in real circumstances, the deterioration of vaccine efficacy is much more complicated. According to a recent study [70], there existed substantial waning immunity after the BNT162b2 vaccine (Pfizer-BioNTech) in Israel while coping with the Delta variant. The research revealed that the vaccine efficacy for Delta variant varied considerably (ranging from 50+% to 80+%) for different age groups that were adjusted for weeks of infection, which verified that the vaccine efficacy for Delta variant could not be simply interpreted by a fixed value. Furthermore, the diffusion of a new variant also requires some time to reach dominance. For the epidemic spreading in Israel, the diffusion of the Delta variant rose from a proportion of less than 10% to nearly 100% in over two months (from April 2021 to June 2021; data source, CoVariants and GISAID), which was demonstrated as an evolving process. Hence, instead of assigning fixed values for pre-Delta and post-Delta vaccine efficacy, advisably, we suggest that mixture models would be more plausible for explaining scenarios of new variants along with different vaccine efficacy, which is surely one of the main topics for our future work.

From a broader perspective, here, we try to explore the following unsolved questions about the mitigation measures. Which containment policy might be the most effective one? Do we have to keep the strictest mitigation strategies until the end of the epidemic? If there emerges a new variant with different transmissibility, will the current containment measures still be effective? Admittedly, it is true that our proposed Bayesian framework model lacks the specific design to sufficiently answer those questions from both the scientific and medical viewpoints; however, some other studies have provided attempts for assessing the efficacy of various containment measures via compartmental models. Sharov [71] presented a susceptible–infected–recovered (SIR)-based model equipped with a system of Kolmogorov linear differential equations. Using epidemiological data as input parameters, the proposed model calculated times and levels of herd immunity formation in 15 countries (13 of which had lockdown measures, while the other two had no-lockdown measures), and the effectiveness of the lockdown measures was evaluated by the comparison of these output parameters. Results showed that only at the initial phase (within 7–8 weeks) of the epidemic, the lockdown mode may bring better results in COVID-19 containment than the no-lockdown one. Universally, the study demonstrated no significant difference between lockdown and no-lockdown modes of COVID-19 containment in terms of times and levels of herd immunity. Valba et al. [72] evaluated the effectiveness of two means of clusterization, namely, the self-isolation and sharp clustering (border closure), by running SIR models on e-networks (evolutionary-grown networks from a randomly generated Erdős–Rényi graph) and i-networks (instantaneously created clustered networks) The study demonstrated that a network, evolutionary grown, blocks the spread of COVID-19 better than an instantly clustered network. When a considerable part of population is already infected, an ideologically loaded mechanism of preventing COVID-19 is more effective than the other. The results proved that mathematics may provide a means of comparing ideologies in public policy, which is crucial for analyzing the future development of societies. Nenchev [73] studied the optimal control of COVID-19 spread based on an augmented SIR model with permanent immunity and no vaccine availability. The proposed model introduced a time varying quarantine strength control $u\left(t\right)$ along with the corresponding quarantined (including the recovered individuals) population Q. Assuming limited isolation control and capacity constraints, an optimal strategy for quarantine strength control that balances between the scale of infections and the effort of isolation was mathematically derived from essential optimality conditions. The numerical simulations indicated that a sufficiently delayed and controlled release of the lockdown are optimal for overcoming the epidemic outbreak.

Above all, these related works have substantially expanded the social analysis, and their methods could technically be adapted to our Bayesian framework for different scenarios. Other improvements include a more explanatory design for the likelihood, additional dynamics of reinfections, and the extension of the third dose of vaccination. Bayesian mixture models could also be supplemented into our framework to depict the different transmissibilities and vaccine efficacies of new variants, as well as the distinctive protection rates of various types of vaccines.

We believe that our modeling framework could be adopted for case studies of other countries which have experienced several waves of epidemic and are highly vaccinated. Additional details, including the varying containment measures implemented by the local government, are critical in consequence of spatial heterogeneity, which may lead to different sets of priors (especially the priors for change points) based on their real situations.