Mathematical Modeling and Short-Term Forecasting of the COVID-19 Epidemic in Bulgaria: SEIRS Model with Vaccination

Abstract

Data from the World Health Organization indicate that Bulgaria has the second-highest COVID-19 mortality rate in the world and the lowest vaccination rate in the European Union. In this context, to find the crucial epidemiological parameters that characterize the ongoing pandemic in Bulgaria, we introduce an extended SEIRS model with time-dependent coefficients. In addition to this, vaccination and vital dynamics are included in the model. We construct an appropriate Cauchy problem for a system of nonlinear ordinary differential equations and prove that its unique solution possesses some biologically reasonable features. Furthermore, we propose a numerical scheme and give an algorithm for the parameters identification in the obtained discrete problem. We show that the found values are close to the parameters values in the original differential problem. Based on the presented analysis, we develop a strategy for short-term prediction of the spread of the pandemic among the host population. The proposed model, as well as the methods and algorithms for parameters identification and forecasting, could be applied to COVID-19 data in every single country in the world.

1. Introduction

Mathematical and computer modelling is of great practical importance in controlling and predicting the spread of various viruses and the development of pandemics caused by them. SEIR-based models are the most used deterministic models for this purpose. In the present paper, we explore a time-dependent generalized SEIRS type model with vaccination and vital dynamics, called the SEIRS-VB model. In this model, the dynamics of the infection in six groups of the host population are modelled by a Cauchy problem for a system of nonlinear ordinary differential equations. The new model is a generalization of the model introduced in our previous work [1].

The pandemic of coronavirus disease COVID-19 caused by the coronavirus SARS-CoV-2 is characterized by significant morbidity and mortality. The novel coronavirus first identified in Wuhan, Hubei province, China at the end of 2019 has spread worldwide and the World Health Organization (WHO) declared five variants of concern (VOC)—three previously circulating VOCs: Alpha, Beta, Gamma, and two currently circulating Delta and Omicron. Several new vaccines were created and being endorsed for emergency use at the end of 2020. The current vaccines may not be a perfect fit for the new Omicron variant, but they are still the best line of defense against COVID-19. The coronavirus pandemic in each country directly affected economic and social life. Consequently, the development of effective mathematical models and methods for predicting the spread and development of the epidemiological process is directly related to many societal challenges.

We note that The National Health Commission of China reported new 3486 symptomatic cases and 20,782 asymptomatic infections with COVID-19 on 15 April 2022. This means that the COVID-19 pandemic is not over yet, and this topic will be relevant in the future as well.

Recently, there has been a lot of interest in making realistic short- and long-term projection for COVID-19 transmission dynamics in a number of countries using deterministic SEIR-based models or creating different possible scenarios for future development of the pandemic depending on the applied vaccination strategies and measures to limit the spread. For example, in [2] a SEIRS model with demography and constant coefficients is used for predicting long-term scenarios and analyzing the time it takes to reach an “endemic equilibrium”. Discrete SIR-type models with constant coefficients are used in [3] to predict the COVID-19 pandemic turning point and ending time in the USA. A prediction method, based on Taylor series expansion for an SIR-like model is applied to Canadian COVID-19 data in [4]. For medium-term and long-term COVID-19 forecasts different optimization algorithms [5] and stochastic methods such as simulated annealing, differential evolution, and genetic algorithm in case of SEIR-D-type models are used (see [6,7,8]). In a number of countries, a change in recovery and latency rates behaviour has been observed during the domination of different variants of SARS-CoV-2 or rapid change of the transmission rate because of government control measures. Hence, it is not realistic to use models with constant parameters over a long period of time. This is why such hypotheses are mostly used for short-term forecasts with SIR models. These predictions are usually made under the assumption that some of the coefficients are constant in a short time interval. In this case, due to the high variability of daily values, different data smoothing algorithms are used. For example, in [9] a short-term prediction methodology is suggested in the case of the classical SIR model, where an infection change ratio is assumed to be a constant for different periods. In [10] Poisson distribution for the daily incidence number, and a gamma distribution for the series interval, are used to estimate the effective reproduction number, which is supposed to stay the same in the forecasting step made by the SIR model. A distributed optimal control epidemiological model is presented in [11] and interesting features of the optimal policy for social distancing are shown. The ensemble Kalman filter as a good short-term predictor is used in [12,13,14], where an algorithm for short-term forecasting based on the estimation of the contact parameter in a stochastic SEIR model with sequential data assimilation is suggested. In [15] we developed a strategy for 14-day prediction of all compartments in the SEIR model, based on the prediction of the parameters daily values in Bulgaria. Our strategy takes into account the level of the government control measures, the new cases/tests dependency on the day of the week, and the length of the healing process. SIR-based models for studying the spread of the COVID-19 pandemic in Bulgaria are also applied in [16,17,18]. Another interesting topic is the modeling of spread and dynamics of COVID-19 with appropriate fractional models. The reason of this is the memory effects of the disease (dependence not only of the current state of infected people, but also from the situation in the past). Let us mention two publications in this area [19,20]. In the first one, a two-side fractional generalized SEIR model is proposed, and the key epidemiological parameters of COVID-19 pandemic in the United States are identified and ranked. In [20] different integer-order and fractional-order models are explored, and their performance with COVID-19 data in China is analyzed.

Recently, several SEIR-based models with vaccination have been studied and used for modeling of COVID-19 pandemic in different countries. In some models [21,22] it is assumed that all vaccinated persons are well protected. While in others [1,3,23,24] as in the present SEIRS-VB model, it is assumed that the vaccine is imperfect and vaccinated persons are not susceptible only for some period of time.

On the other hand, some discrete models are very often used instead of the differential deterministic models (see [3,12,18,25,26]). We will note that the discrete and the differential models can be considered as “equivalent” only if the step size in the discrete model is small enough. Actually, this paper is looking for a discrete model that can be used in the case of real data on the one hand (i.e., with a one-day step for active cases, for example). On the other hand, it is close enough (in appropriate sense) to the new differential model SEIRS-VB. The new model is used for modeling and short-term forecasting the spread of the COVID-19 pandemic in Bulgaria by considering publicly available data from 8 March 2020 to 10 April 2022.

The present paper is organized as follows:

Section 2 begins with the formulation of new SEIRS-VB problem. In Section 2.2, we study the analytic properties of the solution of this problem and show that it has biologically reasonable properties. We prove the existence of a unique solution of the problem, which is well defined and non-negative. In Section 3, we introduce a semi-implicit discrete model and prove that it has similar biological properties like the differential model SEIRS-VB if the step size is small enough. Since available data sets contain daily values of some functions in the differential model, we can assume that the step size in the discrete model is a priori fixed. We divide the parameters into two groups: (i) those that can be selected from the available statistical data and (ii) unknown ones (which should be found using the model). In this context, in Section 4, we formulate an appropriate “inverse” discrete problem IDP to find the unknown values in the discrete model. Furthermore, we present an algorithm for solving this problem. Using available COVID-19 data, in Section 5, the theoretical study is used for identification of the parameters in the discrete model and for experimental analysis of the COVID-19 pandemic situation in Bulgaria. Comparison with the model, in Section 6, using the obtained values of parameters, we solve numerically the original differential problem and show that it gives good approximation of the number of the active cases in norms of the $\left\{l_p\right\}$ family. Based on these results, in Section 7, we propose a forecasting methodology and make numerical experiments for prediction of the numbers of: the active cases, the new daily cases, the cumulative COVID-19 deaths, and the cumulative number of the recovered individuals in Bulgaria for two different 14-day periods. Discussion and Conclusions are outlined in Section 8 and Section 9, respectively.

2. Time-Dependent SEIRS-Based Model with Vaccination and Vital Dynamics

2.1. The SEIRS-VB Model: Formulation

To formulate the new model, named SEIRS-VB, the host population is divided into six compartments, as follows (see Figure 1):

• $S\left(t\right)$—Susceptible individuals. These are the people who may be infected and can become virus carriers. In this group, we include all individuals without immunity (unvaccinated, not fully vaccinated, vaccinated people for whom the vaccine is ineffective, fully vaccinated or recovered individuals who have lost their immunity). Usually at the beginning of a pandemic, as in the case of COVID-19, the whole host population is susceptible.
• $E\left(t\right)$—Exposed individuals. These are virus carrier individuals in the latent stage, during which they are not virus spreaders. They usually have no symptoms.
• $I\left(t\right)$—Infectious individuals. These are virus carriers and virus spreaders of extremely high infectivity. The former are likely to transmit the virus in case of contact.
• $R\left(t\right)$—Recovered individuals with immunity. These individuals have disease acquired immunity. They have recovered, and thus are protected from the disease.
• $V\left(t\right)$—Vaccinated susceptible individuals. These are fully vaccinated persons for whom the vaccine is effective. However, they have not developed antibodies. They can do so after a certain period of time or else they will become exposed individuals before that. It is worth pointing out that, due to the vaccine imperfection, some of the vaccinated individuals can not develop antibodies, and they can not pass from group $S\left(t\right)$ to group $V\left(t\right)$.
• $B\left(t\right)$—Individuals with vaccination-acquired immunity. These are vaccinated individuals who are well protected from future infection because they have antibodies.

Figure 1. The model diagram.

Figure 1. The model diagram.

SEIRS-VB model (Figure 1) involves three “chains” S-E-I-R-S, S-V-B-S and S-V-E-I-R-S.

(1)

The S-E-I-R-S chain (horizontal red line in Figure 1) describes infection transmission among unvaccinated individuals. Each susceptible individual from (S) can be infected in case of contact with infectious persons and is transferred to group (E) with transmission rate $\beta \left(t\right).$ Furthermore, the individual goes in group (I) with latency rate $\omega \left(t\right)$ and later to group (R) with recovery rate $\gamma \left(t\right)$. Later, this person loses disease-acquired immunity and moves again to group (S) with re-infection rate $\lambda \left(t\right)$.

(2)

The chains S-V-B-S and S-V-E-I-R-S describe the change in the proportion of the numbers of the groups that contain vaccinated persons. Each vaccinated person in (S) can move to group (V) with vaccination parameter $\alpha \left(t\right)$ or stay in group (S) due to imperfect vaccines which means that this does not provide 100 % safety (with effectiveness $\sigma <1$). Each individual in group (V) can develop antibodies and go to group (B) with antibodies rate $\mu \left(t\right)$ or be infected before that and move to group (E) with the transmission rate $\beta \left(t\right)$. Individuals in group (B) lose vaccination-acquired immunity and move again to group (S) with re-infection rate $\nu \left(t\right)$.

(3)

The SEIRS-VB model describes vital dynamics as well. We assume that all new born individuals are susceptible and they come to group (S) with birth rate $\Lambda \left(t\right)$. At the same time, individuals in each model’s group can move out from the model, due to natural mortality with rate $\theta \left(t\right)$. Of course, infectious individuals (I) can leave the model due to COVID-19 mortality with rate $\tau \left(t\right)$. Since we use a model with re-susceptibility, it is reasonable to assume that all individuals having received booster doses of vaccines are new fully vaccinated persons, i.e., belong to group (S) who can possibly move to (V) due to these doses.

It is natural that the total population size is the number of all living individuals

$$N\left(t\right):=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)+B\left(t\right)+V\left(t\right).$$

The SEIRS-VB model is described by the following Cauchy problem for a system of nonlinear ordinary differential equations

$$\left|\begin{array}{c}{\displaystyle \frac{dS}{dt}=\Lambda \left(t\right)N\left(t\right)-(\alpha \left(t\right)+\theta \left(t\right))S\left(t\right)-\frac{\beta \left(t\right)}{N\left(t\right)}S\left(t\right)I\left(t\right)+\lambda \left(t\right)R\left(t\right)+\nu \left(t\right)B\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dE}{dt}=\frac{\beta \left(t\right)}{N\left(t\right)}S\left(t\right)I\left(t\right)+\frac{\beta \left(t\right)}{N\left(t\right)}V\left(t\right)I\left(t\right)-(\omega \left(t\right)+\theta \left(t\right))E\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dI}{dt}=\omega \left(t\right)E\left(t\right)-(\gamma \left(t\right)+\tau \left(t\right)+\theta \left(t\right))I\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dR}{dt}=\gamma \left(t\right)I\left(t\right)-(\lambda \left(t\right)+\theta \left(t\right))R\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dV}{dt}=\alpha \left(t\right)S\left(t\right)-(\mu \left(t\right)+\theta \left(t\right))V\left(t\right)-\frac{\beta \left(t\right)}{N\left(t\right)}I\left(t\right)V\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dB}{dt}=\mu \left(t\right)V\left(t\right)-(\nu \left(t\right)+\theta \left(t\right))B\left(t\right)}\hfill \end{array}\right.$$

(1)

with non-negative initial conditions

$$S\left(t_0\right)=S_0,E\left(t_0\right)=E_0,I\left(t_0\right)=I_0,R\left(t_0\right)=R_0,V\left(t_0\right)=V_0,B\left(t_0\right)=B_0,$$

(2)

where $t_0\ge 0$ is a real number.

The coefficients in the system (1) are time-dependent, and they involve the parameters listed in

Table 1. Parameters of the SEIRS-VB model.

Parameter	Description	Units
$\Lambda \left(t\right)$	birth rate	$\frac{births}{population}$$/day$
$a\left(t\right)$	vaccination rate	$\frac{vaccinated}{population}/day$
$\sigma$	vaccine effectiveness	$\frac{excessrisk}{riskamongvaccinated}$
$\alpha \left(t\right)$	vaccination parameter	$\sigma a\left(t\right)$
$\beta \left(t\right)$	transmission rate	1/days
$\gamma \left(t\right)$	recovery rate	1/ days
$\omega \left(t\right)$	latency rate	1/days
$\theta \left(t\right)$	natural mortality rate	$\frac{deaths}{population}/day$
$\tau \left(t\right)$	mortality rate of infectious people	$\frac{deaths}{infectious}/day$
$\lambda \left(t\right)$	reinfection rate of recovered individuals	1/days
$\nu \left(t\right)$	reinfection rate of vaccinated individuals	1/days
$\mu \left(t\right)$	antibody rate	1/days

.

In the classical SIR and SEIR models [27], as well as in many SIR/SEIR-based models with constant coefficients, it is assumed that the population size N is a constant. In the model (1), the population size $N\left(t\right)$ is a time-dependent function. Summing up the equations in the system (1), we obtain

$$\frac{dN}{dt}=[\Lambda \left(t\right)-\theta \left(t\right)]N\left(t\right)-\tau \left(t\right)I\left(t\right).$$

(3)

2.2. Properties of the Analytic Solution of the SEIRS-VB Model

To study the biologically reasonable futures of model in the the case when the population is fixed (1), following [28,29,30] we introduce the functions

$$\begin{array}{cc}\hfill & \tilde{S}\left(t\right):=\frac{S\left(t\right)}{N\left(t\right)},\tilde{E}\left(t\right):=\frac{E\left(t\right)}{N\left(t\right)},\tilde{I}\left(t\right):=\frac{I\left(t\right)}{N\left(t\right)},\hfill \\ \hfill & \tilde{R}\left(t\right):=\frac{R\left(t\right)}{N\left(t\right)},\tilde{V}\left(t\right):=\frac{V\left(t\right)}{N\left(t\right)},\tilde{B}\left(t\right):=\frac{B\left(t\right)}{N\left(t\right)},\hfill \end{array}$$

(4)

where, obviously, $N\left(t\right)>0$ for the considered time-frame $[t_0,T]$.

For the new unknowns and from (1) and (2), we obtain the following reduced Cauchy problem:

$$\left|\begin{array}{c}{\displaystyle \frac{d\tilde{S}}{dt}=\Lambda \left(t\right)-(\alpha \left(t\right)+\Lambda \left(t\right))\tilde{S}\left(t\right)+(\tau \left(t\right)-\beta \left(t\right))\tilde{S}\left(t\right)\tilde{I}\left(t\right)+\lambda \left(t\right)\tilde{R}\left(t\right)+\nu \left(t\right)\tilde{B}\left(t\right),}\hfill \\ \\ {\displaystyle \frac{d\tilde{E}}{dt}=\beta \left(t\right)\tilde{S}\left(t\right)\tilde{I}\left(t\right)+\beta \left(t\right)\tilde{V}\left(t\right)\tilde{I}\left(t\right)-(\omega \left(t\right)+\Lambda \left(t\right))\tilde{E}\left(t\right)+\tau \left(t\right)\tilde{E}\left(t\right)\tilde{I}\left(t\right),}\hfill \\ \\ {\displaystyle \frac{d\tilde{I}}{dt}=\omega \left(t\right)\tilde{E}\left(t\right)-(\gamma \left(t\right)+\tau \left(t\right)+\Lambda \left(t\right))\tilde{I}\left(t\right)+\tau \left(t\right){\left(\tilde{I}\left(t\right)\right)}^2,}\hfill \\ \\ {\displaystyle \frac{d\tilde{R}}{dt}=\gamma \left(t\right)\tilde{I}\left(t\right)-(\lambda \left(t\right)+\Lambda \left(t\right))\tilde{R}\left(t\right)+\tau \left(t\right)\tilde{R}\left(t\right)\tilde{I}\left(t\right),}\hfill \\ \\ {\displaystyle \frac{d\tilde{V}}{dt}=\alpha \left(t\right)\tilde{S}\left(t\right)-(\mu \left(t\right)+\Lambda \left(t\right))\tilde{V}\left(t\right)+(\tau \left(t\right)-\beta \left(t\right))\tilde{I}\left(t\right)\tilde{V}\left(t\right),}\hfill \\ \\ {\displaystyle \frac{d\tilde{B}}{dt}=\mu \left(t\right)\tilde{V}\left(t\right)-(\nu \left(t\right)+\Lambda \left(t\right))\tilde{B}\left(t\right)+\tau \left(t\right)\tilde{B}\left(t\right)\tilde{I}\left(t\right)}\hfill \end{array}\right.$$

(5)

with non-negative initial conditions

$$\tilde{S}\left(t_0\right)={\tilde{S}}_0,\tilde{E}\left(t_0\right)={\tilde{E}}_0,\tilde{I}\left(t_0\right)={\tilde{I}}_0,\tilde{R}\left(t_0\right)={\tilde{R}}_0,\tilde{V}\left(t_0\right)={\tilde{V}}_0,\tilde{B}\left(t_0\right)={\tilde{B}}_0,$$

(6)

where

$${\tilde{S}}_0+{\tilde{E}}_0+{\tilde{I}}_0+{\tilde{R}}_0+{\tilde{V}}_0+{\tilde{B}}_0=1.$$

(7)

Summing up the equations in (5), for

$$\tilde{N}\left(t\right):=\tilde{S}\left(t\right)+\tilde{E}\left(t\right)+\tilde{I}\left(t\right)+\tilde{R}\left(t\right)+\tilde{B}\left(t\right)+\tilde{V}\left(t\right)$$

we obtain

$$\frac{d\tilde{N}}{dt}=0$$

and because of (7) obviously

$$\tilde{N}\left(t\right)\equiv 1.$$

Let us introduce the notations

$$\tilde{x}\left(t\right):=\left(\tilde{S}\left(t\right),\tilde{E}\left(t\right),\tilde{I}\left(t\right),\tilde{R}\left(t\right),\tilde{V}\left(t\right),\tilde{B}\left(t\right)\right),$$

$${\tilde{x}}_0:=\left({\tilde{S}}_0,{\tilde{E}}_0,{\tilde{I}}_0,{\tilde{R}}_0,{\tilde{V}}_0,{\tilde{B}}_0\right),$$

$$p\left(t\right):=(\Lambda (t),\theta (t),\omega (t),\mu (t),\alpha (t),\beta (t),\gamma (t),\tau (t),\nu (t),\lambda (t\left)\right)$$

(8)

and

$$f(t,\tilde{x})=\left(\begin{array}{c}\Lambda \left(t\right)-(\alpha \left(t\right)+\Lambda \left(t\right))\tilde{S}+(\tau \left(t\right)-\beta \left(t\right))\tilde{S}\tilde{I}+\lambda \left(t\right)\tilde{R}+\nu \left(t\right)\tilde{B}\\ \beta \left(t\right)\tilde{S}\tilde{I}+\beta \left(t\right)\tilde{V}\tilde{I}-(\omega \left(t\right)+\Lambda \left(t\right))\tilde{E}+\tau \left(t\right)\tilde{E}\tilde{I}\\ \omega \left(t\right)\tilde{E}-(\gamma \left(t\right)+\tau \left(t\right)+\Lambda \left(t\right))\tilde{I}+\tau \left(t\right){\tilde{I}}^2\\ \gamma \left(t\right)\tilde{I}-(\lambda \left(t\right)+\Lambda \left(t\right))\tilde{R}+\tau \left(t\right)\tilde{R}\tilde{I}\\ \alpha \left(t\right)\tilde{S}-(\mu \left(t\right)+\Lambda \left(t\right))\tilde{V}-\beta \left(t\right)\tilde{V}\tilde{I}+\tau \left(t\right)\tilde{V}\tilde{I}\\ \mu \left(t\right)\tilde{V}-(\nu \left(t\right)+\Lambda \left(t\right))\tilde{B}+\tau \left(t\right)\tilde{B}\tilde{I}\end{array}\right).$$

(9)

Then, the Cauchy problem (5) and (6) can be rewritten in the form

$$\left|\begin{array}{c}\dot{\tilde{x}}\left(t\right)=f(t,\tilde{x}),\hfill \\ \tilde{x}\left(t_0\right)={\tilde{x}}_0.\hfill \end{array}\right.$$

(10)

Remark 1.

According to the notation: here and further, we write that a vector is nonnegative/positive if all its components are nonnegative/positive.

Remark 2.

Let us note that the initial conditions for the pandemic diffusion models with $N=1$ usually satisfy the follows assumptions (see also [28]):

$$0<{\tilde{I}}_0\ll 1,0\le {\tilde{V}}_0\ll 1,0\le {\tilde{E}}_0\ll 1,{\tilde{S}}_0=1-{\tilde{E}}_0-{\tilde{I}}_0-{\tilde{V}}_0,{\tilde{R}}_0=0,{\tilde{B}}_0=0.$$

Therefore, the conditions of the following theorem are natural.

Theorem 1.

Let ${\tilde{x}}_0\ge 0$, ${\tilde{S}}_0>0,{\tilde{I}}_0>0$, $p\in C\left([t_0,T]\right)$ and $p\left(t\right)\ge 0,\Lambda \left(t\right)>0,\beta \left(t\right)>0,\omega \left(t\right)>0$ for $t_0\le t\le T$. Then, there exists a unique solution $\tilde{x}\left(t\right)$ of the Cauchy problem (10), which is well defined and bounded for all $t\in [t_0,T]$ and $\tilde{x}\left(t\right)\ge 0,$$\tilde{S}\left(t\right)>0,\tilde{I}\left(t\right)>0$ for $t\in [t_0,T].$

Proof. 

According to (8), the nonnegativity of $p\left(t\right)$ in the assumption of this theorem means that $\Lambda \left(t\right)\ge 0,\theta \left(t\right)\ge 0,\omega \left(t\right)\ge 0,\mu \left(t\right)\ge 0,\alpha \left(t\right)\ge 0,\beta \left(t\right)\ge 0,\gamma \left(t\right)\ge 0,\tau \left(t\right)\ge 0,\nu \left(t\right)\ge 0,\lambda \left(t\right)\ge 0$. In view of Table 1, the nonnegativity of the parameters in $p\left(t\right)$ is a biologically reasonable property. Let us note first that, from the special kind of the function $f(t,\tilde{x})$ (second-order polynomials with respect to $\tilde{x}$, see (9)), it follows that $f,\frac{\partial f}{\partial {\tilde{x}}_j}\in C([t_0,T]\times {\mathbb{R}}^6),j=1,2,\cdots ,6$. According to Picard’s existence and uniqueness theorem (Theorem 1.1., p. 8 in [31]), there exists a unique solution $\tilde{x}\left(t\right)$ of the Cauchy problem (10) defined in the interval $[t_0,T_1]$ for some $T_1\in (t_0,T]$. Now, we would like to show that $T_1=T$.

Now, we prove that $\tilde{x}\left(t\right)\ge 0$ and $\tilde{S}\left(t\right)>0,\tilde{I}\left(t\right)>0$ for $t\in [t_0,T_1]$. First assume, on the contrary, that $\tilde{S}\left(t\right)\tilde{I}\left(t\right)$ vanishes in $(t_0,T_1]$. Denote with $t_1\in (t_0,T_1]$ the first time such that $\tilde{S}\left(t_1\right)\tilde{I}\left(t_1\right)=0,$ i.e., $\tilde{S}\left(t\right)>0$ and $\tilde{I}\left(t\right)>0$ for $t\in [t_0,t_1)$.

Then, for $t\in [t_0,t_1]$, we have from the system (5):

$$\tilde{R}\left(t\right)=e^{{\int }_{t_0}^t[\tau \left(s\right)\tilde{I}\left(s\right)-\lambda \left(s\right)-\Lambda \left(s\right)]ds}\left({\tilde{R}}_0+{\int }_{t_0}^t\gamma \left(s\right)\tilde{I}\left(s\right)e^{{\int }_{t_0}^s[\lambda \left(\sigma \right)+\Lambda \left(\sigma \right)-\tau \left(\sigma \right)\tilde{I}\left(\sigma \right)]d\sigma }ds\right)\ge 0,$$

(11)

$$\begin{array}{cc}\hfill \tilde{V}\left(t\right)=& e^{-{\int }_{t_0}^t[\mu \left(s\right)+\Lambda \left(s\right)+(\beta \left(s\right)-\tau \left(s\right))\tilde{I}\left(s\right)]ds}\hfill \\ \hfill & \times \left({\tilde{V}}_0+{\int }_{t_0}^t\alpha \left(s\right)\tilde{S}\left(s\right)e^{{\int }_{t_0}^s[\mu \left(\sigma \right)+\Lambda \left(\sigma \right)+(\beta \left(\sigma \right)-\tau \left(\sigma \right))\tilde{I}\left(\sigma \right)]d\sigma }ds\right)\ge 0,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \tilde{E}\left(t\right)=& e^{-{\int }_{t_0}^t[\omega \left(s\right)+\Lambda \left(s\right)-\tau \left(s\right)\tilde{I}\left(s\right)]ds}\hfill \\ \hfill & \times \left({\tilde{E}}_0+{\int }_{t_0}^t\beta \left(s\right)[\tilde{S}\left(s\right)+\tilde{V}\left(s\right)]\tilde{I}\left(s\right)e^{{\int }_{t_0}^s[\omega \left(\sigma \right)+\Lambda \left(\sigma \right)-\tau \left(\sigma \right)\tilde{I}\left(\sigma \right)]d\sigma }ds\right)>0,\hfill \end{array}$$

(13)

$$\tilde{B}\left(t\right)=e^{{\int }_{t_0}^t[\tau \left(s\right)\tilde{I}\left(s\right)-\nu \left(s\right)-\Lambda \left(s\right)]ds}\left({\tilde{B}}_0+{\int }_{t_0}^t\mu \left(s\right)\tilde{V}\left(s\right)e^{{\int }_{t_0}^s[\nu \left(\sigma \right)+\Lambda \left(\sigma \right)-\tau \left(\sigma \right)\tilde{I}\left(\sigma \right)]d\sigma }ds\right)\ge 0.$$

(14)

Note that, for the inequalities (13) and (14) we had used the already proved (11) and (12).

Now, from the first equation of the system (5), we have

$$\begin{array}{cc}\hfill \tilde{S}\left(t_1\right)=& e^{{\int }_{t_0}^{t_1}[(\tau \left(s\right)-\beta \left(s\right))\tilde{I}\left(s\right)-\alpha \left(s\right)-\Lambda \left(s\right)]ds}\hfill \\ \hfill & \times \left({\tilde{S}}_0+{\int }_{t_0}^{t_1}[\Lambda \left(s\right)+\lambda \left(s\right)\tilde{R}\left(s\right)+\nu \left(s\right)\tilde{B}\left(s\right)]e^{{\int }_{t_0}^s[(\beta \left(\sigma \right)-\tau \left(\sigma \right))\tilde{I}\left(\sigma \right)+\alpha \left(\sigma \right)+\Lambda \left(\sigma \right)]d\sigma }ds\right)>0,\hfill \end{array}$$

because of (12)–(14) and ${\tilde{S}}_0>0$.

Therefore, $\tilde{I}\left(t_1\right)=0$. Then, from the third equation of the system (5), we have

$$\frac{d\tilde{I}}{dt}\left(t_1\right)=\omega \left(t_1\right)\tilde{E}\left(t_1\right)>0,$$

and, because $\tilde{I}\left(t\right)>0$ for $t_0\le t<t_1$, it leads to a contradiction with $\tilde{I}\left(t_1\right)=0$. It follows that $\tilde{I}\left(t\right)$ is also always positive for $t\in [t_0,T_1]$.

Hence, $\tilde{S}\left(t\right)>0$ and $\tilde{I}\left(t\right)>0$ for $t\in [t_0,T_1]$. Now, it is easy to see that the inequalities (11)–(14) hold for all $t\in [t_0,T_1]$. In this way, we show that $\tilde{x}\left(t\right)\ge 0$ for $t\in [t_0,T_1]$.

Since $\tilde{x}\left(t\right)\ge 0$ and $\tilde{S}\left(t\right)+\tilde{E}\left(t\right)+\tilde{I}\left(t\right)+\tilde{R}\left(t\right)+\tilde{V}\left(t\right)+\tilde{B}\left(t\right)=1$, we obtain that $\parallel \tilde{x}\parallel \le \sqrt{6}$, i.e., the solution is uniformly bounded on $[t_0,T_1]$. Now, by Extension theorem (Theorem 3.1., p. 12 in [31]), we conclude that $T_1\equiv T,$ i.e., the solution $\tilde{x}\left(t\right)$ exists on $[t_0,T]$. This completes the proof. □

To study the solutions of original problem (1) and (2), we introduce the notations

$$x\left(t\right):=\left(S\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right),V\left(t\right),B\left(t\right)\right),$$

$$x_0:=\left(S_0,E_0,I_0,R_0,V_0,B_0\right).$$

Theorem 2.

Let $x_0\ge 0$, $S_0>0,I_0>0$, $p\in C\left([t_0,T]\right)$ and $p\left(t\right)\ge 0,\Lambda \left(t\right)>0,\beta \left(t\right)>0,\omega \left(t\right)>0$ for $t_0\le t\le T$. Then, there exists a solution $x\left(t\right)$ of the Cauchy problem (1) and (2), which is defined for $t\in [t_0,T]$ and $x\left(t\right)\ge 0,$$S\left(t\right)>0,I\left(t\right)>0$ for $t\in [t_0,T].$ The solution with such properties is unique.

Proof. 

To prove the existence and uniqueness of such a solution, we use Theorem 1.

(i) Existence. Denote

$$N_0:=S_0+E_0+I_0+R_0+V_0+B_0>0$$

(because $x_0\ge 0$ with $S_0>0$ and $I_0>0$) and let us consider the Cauchy problem (10) with initial condition ${\tilde{x}}_0:=x_0/N_0\ge 0.$

Since ${\tilde{S}}_0=S_0/N_0>0$, ${\tilde{I}}_0=I_0/N_0>0$ and all other conditions of Theorem 1 are fulfilled, Theorem 1 gives a function $\tilde{x}\left(t\right)=(\tilde{S}\left(t\right),$$\tilde{E}\left(t\right),\tilde{I}\left(t\right),\tilde{R}\left(t\right),\tilde{V}\left(t\right),\tilde{B}\left(t\right))\ge 0$, which is a unique solution of the problem (10), defined for $t\in [t_0,T]$.

Let us consider the Cauchy problem

$$\left|\begin{array}{c}\dot{N}\left(t\right)=[\Lambda \left(t\right)-\theta \left(t\right)-\tau \left(t\right)\tilde{I}\left(t\right)]N\left(t\right),\hfill \\ N\left(t_0\right)=N_0.\hfill \end{array}\right.$$

(15)

Obviously, the unique solution of problem (15) is

$$N\left(t\right)=N_0{e}^{{\int }_{t_0}^t[\Lambda \left(s\right)-\theta \left(s\right)-\tau \left(s\right)\tilde{I}\left(s\right)]ds}>0,$$

and it is defined for $t\in [t_0,T]$.

Above, we transformed the system (1) to the system (5). Now, following the reverse way and using $N\left(t\right)>0$, we obtain that the function $x\left(t\right):=N\left(t\right)\tilde{x}\left(t\right)\ge 0$ is a solution of the Cauchy problem (1) and (2) in $[t_0,T]$ and $S\left(t\right)>0,I\left(t\right)>0$ for all $t\in [t_0,T]$.

(ii) Uniqueness. Let

$$x^1\left(t\right)=(S^1\left(t\right),E^1\left(t\right),I^1\left(t\right),R^1\left(t\right),B^1\left(t\right),V^1\left(t\right))$$

and

$$x^2\left(t\right)=(S^2\left(t\right),E^2\left(t\right),I^2\left(t\right),R^2\left(t\right),B^2\left(t\right),V^2\left(t\right))$$

be two solutions of the problem (1) and (2) with initial condition $x^1\left(t_0\right)=x^2\left(t_0\right)=x_0$, which are defined for $t\in [t_0,T]$ and $x^1\left(t\right)\ge 0,x^2\left(t\right)\ge 0,$$S^1\left(t\right)>0,I^1\left(t\right)>0,$$S^2\left(t\right)>0,I^2\left(t\right)>0$ for all $t\in [t_0,T]$. We define the functions

$$\begin{array}{cc}\hfill N^1\left(t\right)& :=S^1\left(t\right)+E^1\left(t\right)+I^1\left(t\right)+R^1\left(t\right)+B^1\left(t\right)+V^1\left(t\right)>0,\hfill \\ \hfill N^2\left(t\right)& :=S^2\left(t\right)+E^2\left(t\right)+I^2\left(t\right)+R^2\left(t\right)+B^2\left(t\right)+V^2\left(t\right)>0\hfill \end{array}$$

for $t\in [t_0,T]$ and obviously $N^1\left(t_0\right)=N^2\left(t_0\right)=:N_0>0$.

Then, the functions ${\tilde{x}}^1\left(t\right):=x^1\left(t\right)/N^1\left(t\right)\ge 0$ and ${\tilde{x}}^2\left(t\right):=x^2\left(t\right)/N^2\left(t\right)\ge 0$ are two solutions of the system (5), defined for $t\in [t_0,T]$ and ${\tilde{x}}^1\left(t_0\right)={\tilde{x}}^2\left(t_0\right)={\tilde{x}}_0:=x_0/N_0\ge 0$ with ${\tilde{S}}^0:=S_0/N_0>0,{\tilde{I}}^0:=I_0/N_0>0$ and ${\tilde{S}}^1\left(t\right)>0,{\tilde{I}}^1\left(t\right)>0,$ ${\tilde{S}}^2\left(t\right)>0,{\tilde{I}}^2\left(t\right)>0$ for $t\in [t_0,T]$. Now, Theorem 1 gives ${\tilde{x}}^1\left(t\right)={\tilde{x}}^2\left(t\right)$ for $t\in [t_0,T]$ and therefore

$$\frac{x^1\left(t\right)}{N^1\left(t\right)}=\frac{x^2\left(t\right)}{N^2\left(t\right)},t_0\le t\le T.$$

(16)

In particular, we have

$$\frac{I^1\left(t\right)}{N^1\left(t\right)}=\frac{I^2\left(t\right)}{N^2\left(t\right)},t_0\le t\le T.$$

Now, summing up the equations in the corresponding systems for $x^1\left(t\right)$ and $x^2\left(t\right)$, respectively, it follows that

$$\begin{array}{cc}\hfill {\dot{N}}^1\left(t\right)& =[\Lambda \left(t\right)-\theta \left(t\right)]N^1\left(t\right)-\tau \left(t\right)I^1\left(t\right),\hfill \\ \hfill {\dot{N}}^2\left(t\right)& =[\Lambda \left(t\right)-\theta \left(t\right)]N^2\left(t\right)-\tau \left(t\right)I^2\left(t\right).\hfill \end{array}$$

Hence,

$$\frac{{\dot{N}}^1\left(t\right)}{N^1\left(t\right)}=\frac{{\dot{N}}^2\left(t\right)}{N^2\left(t\right)},$$

which implies $N^1\left(t\right)/N^2\left(t\right)=const.$ for $t\in [t_0,T].$ Finally, since $N^1\left(t_0\right)=N^2\left(t_0\right)=N_0,$ we conclude that $N^1\left(t\right)=N^2\left(t\right)$ for $t\in [t_0,T]$. Now, (16) gives $x^1\left(t\right)=x^2\left(t\right)$ for $t\in [t_0,T]$.

This completes the proof. □

3. Discretization of the SEIRS-VB Model

In this section, we introduce a discrete analogue of the differential SEIRS-VB model. We are going to show that, if the step size in the discrete model is sufficiently small, the basic properties of its solution are close to those of the solution of the differential model.

Let us denote by

$$A\left(t\right)=I\left(t\right)+E\left(t\right)$$

the number of active cases at the time t, i.e., the number of people who are virus carriers. Some of them are asymptomatic (E), but others are infectious with symptoms (I). The PCR tests can detect RNA from SARS-CoV-2, roughly 1–3 days before the onset of the symptoms. That varies among different variants of virus: 5-14 days for the Wuhan variant, 5–10 for the Alpha variant and 5–7 for the Delta variant and Omicron variant. This means that the virus can be detected no earlier than 48–72 h after the exposure [32]. We assume that the active cases (A) are positive tested individuals, reported by the World Health Organization.

Summing up the second and the third equations in (1), we obtain the following form of SEIRS-VB model:

$$\left|\begin{array}{c}{\displaystyle \frac{dS}{dt}=\Lambda \left(t\right)N\left(t\right)-[\alpha \left(t\right)+\theta \left(t\right)]S\left(t\right)-\frac{\beta \left(t\right)}{N\left(t\right)}S\left(t\right)I\left(t\right)+\lambda \left(t\right)R\left(t\right)+\nu \left(t\right)B\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dA}{dt}=\frac{\beta \left(t\right)}{N\left(t\right)}[S\left(t\right)+V\left(t\right)]I\left(t\right)-\theta \left(t\right)A\left(t\right)-[\gamma \left(t\right)+\tau \left(t\right)]I\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dI}{dt}=\omega \left(t\right)[A\left(t\right)-I\left(t\right)]-[\gamma \left(t\right)+\tau \left(t\right)+\theta \left(t\right)]I\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dR}{dt}=\gamma \left(t\right)I\left(t\right)-[\lambda \left(t\right)+\theta \left(t\right)]R\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dV}{dt}=\alpha \left(t\right)S\left(t\right)-\frac{\beta \left(t\right)}{N\left(t\right)}I\left(t\right)V\left(t\right)-[\mu \left(t\right)+\theta \left(t\right)]V\left(t\right),}\hfill \\ \\ {\displaystyle \frac{dB}{dt}=\mu \left(t\right)V\left(t\right)-[\nu \left(t\right)+\theta \left(t\right)]B\left(t\right)}\hfill \end{array}\right.$$

(17)

with non-negative initial data

$$S\left(t_0\right)=S_0,A\left(t_0\right)=A_0,I\left(t_0\right)=I_0,R\left(t_0\right)=R_0,V\left(t_0\right)=V_0,B\left(t_0\right)=B_0.$$

(18)

Therefore, the total population size is

$$N\left(t\right)=S\left(t\right)+A\left(t\right)+R\left(t\right)+V\left(t\right)+B\left(t\right).$$

We consider the time-frame $t_1,t_2,\cdots ,t_K$, where $t_0\le t_1<t_2<...<t_K\le T$ and introduce the notation for the values of functions

$$(S_k,A_k,I_k,R_k,V_k,B_k)=(S\left(t_k\right),A\left(t_k\right),I\left(t_k\right),R\left(t_k\right),V\left(t_k\right),B\left(t_k\right)),$$

and

$$N_k:=N\left(t_k\right)=S_k+A_k+R_k+V_k+B_k,$$

for $k=1,2,\cdots ,K$ and the values of the parameters

$$\begin{array}{cc}\hfill p_k& :=({\Lambda }_k,{\theta }_k,{\omega }_k,{\mu }_k,{\alpha }_k,{\beta }_k,{\gamma }_k,{\tau }_k,{\nu }_k,{\lambda }_k)\hfill \\ & =(\Lambda \left(t_k\right),\theta \left(t_k\right),\omega \left(t_k\right),\mu \left(t_k\right),\alpha \left(t_k\right),\beta \left(t_k\right),\gamma \left(t_k\right),\tau \left(t_k\right),\nu \left(t_k\right),\lambda \left(t_k\right)),\hfill \end{array}$$

$k=1,2,\cdots ,K-1.$

We introduce the following semi-implicit finite difference scheme as discretization of the Cauchy problem for system (17) considered for $t_1\le t\le T$ with given initial data $S\left(t_1\right)=S_1,A\left(t_1\right)=A_1,I\left(t_1\right)=I_1,R\left(t_1\right)=R_1,V\left(t_1\right)=V_1,B\left(t_1\right)=B_1$:

$$\left|\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \frac{S_k-S_{k-1}}{t_k-t_{k-1}}=}& {\Lambda }_{k-1}{N}_{k-1}-[{\alpha }_{k-1}+{\theta }_{k-1}]S_{k-1}-\frac{{\beta }_{k-1}}{N_{k-1}}{I}_k{S}_{k-1}+{\lambda }_{k-1}{R}_{k-1}+{\nu }_{k-1}{B}_{k-1},\hfill \\ \\ \hfill {\displaystyle \frac{A_k-A_{k-1}}{t_k-t_{k-1}}=}& \frac{{\beta }_{k-1}}{N_{k-1}}[S_{k-1}+V_{k-1}]I_k-{\theta }_{k-1}{A}_{k-1}-[{\gamma }_{k-1}+{\tau }_{k-1}]I_{k-1},\hfill \\ \\ \hfill {\displaystyle \frac{I_k-I_{k-1}}{t_k-t_{k-1}}=}& {\omega }_{k-1}(A_{k-1}-I_{k-1})-[{\gamma }_{k-1}+{\tau }_{k-1}+{\theta }_{k-1}]I_{k-1},\hfill \\ \\ \hfill {\displaystyle \frac{R_k-R_{k-1}}{t_k-t_{k-1}}=}& {\gamma }_{k-1}{I}_{k-1}-[{\lambda }_{k-1}+{\theta }_{k-1}]R_{k-1},\hfill \\ \\ \hfill {\displaystyle \frac{V_k-V_{k-1}}{t_k-t_{k-1}}=}& {\alpha }_{k-1}{S}_{k-1}-\frac{{\beta }_{k-1}}{N_{k-1}}{I}_k{V}_{k-1}-[{\mu }_{k-1}+{\theta }_{k-1}]V_{k-1},\hfill \\ \\ \hfill {\displaystyle \frac{B_k-B_{k-1}}{t_k-t_{k-1}}=}& {\mu }_{k-1}{V}_{k-1}-[{\nu }_{k-1}+{\theta }_{k-1}]B_{k-1},\hfill \end{array}\hfill \end{array}\right.$$

(19)

where $k=2,3,\dots ,K$.

Summing up the following equations: S-equation (the first), A-equation (the second), R-equation (the fourth), V-equation (the fifth) and B-equation (the sixth) in system (19), we obtain

$$\frac{N_k-N_{k-1}}{t_k-t_{k-1}}=[{\Lambda }_{k-1}-{\theta }_{k-1}]N_{k-1}-{\tau }_{k-1}{I}_{k-1}.$$

(20)

Equation (20) can be considered as a discretization of differential Equation (3).

Our first aim is to show that the discrete problem (19) with appropriate initial data has biologically reasonable features, similar to those of the differential problem (1) and (2).

Theorem 3.

Let $S_1>0,A_1\ge I_1>0,R_1\ge 0,V_1\ge 0,B_1\ge 0$, $p_{k-1}\ge 0,{\Lambda }_{k-1}>0,{\beta }_{k-1}>0,{\omega }_{k-1}>0$ for $k=2,3,\cdots ,K$. Let also $t_k-t_{k-1}=h>0$ for all $k=2,3,\cdots ,K,$ and

$$\underset{k=2,\cdots ,K}{max}{q}_{k-1}\le 1/h,$$

(21)

where

$$q_{k-1}:={\theta }_{k-1}+max\{{\alpha }_{k-1}+{\beta }_{k-1}-{\Lambda }_{k-1},{\mu }_{k-1}+{\beta }_{k-1},{\gamma }_{k-1}+{\tau }_{k-1}+{\omega }_{k-1},{\lambda }_{k-1},{\nu }_{k-1}\}.$$

Then, for the values calculated by (19), we have $S_k>0,A_k\ge I_k>0,R_k\ge 0,V_k\ge 0,B_k\ge 0$ for all $k=1,2,\cdots ,K$.

Proof. 

The statement of the theorem holds for $k=1$, and we suppose that $S_{k-1}>0,R_{k-1}\ge 0,V_{k-1}\ge 0,B_{k-1}\ge 0,$$A_{k-1}\ge I_{k-1}>0$ for some k. Therefore, $N_{k-1}>S_{k-1}>0.$

Now, we solve the equations in (19) with respect to $S_k,A_k,I_k,R_k,B_k$ and $V_k$, respectively. Then, since $t_k-t_{k-1}=h$ for all $k=2,3,\cdots ,K,$ and (21) holds, from I-equation (the third), R-equation (the fourth) and B-equation (the sixth) in system (19), it follows that

$$\begin{array}{cc}\hfill I_k& =h{\omega }_{k-1}{A}_{k-1}+[1-h({\gamma }_{k-1}+{\omega }_{k-1}+{\tau }_{k-1}+{\theta }_{k-1})]I_{k-1}>0,\hfill \\ \hfill R_k& =h{\gamma }_{k-1}{I}_{k-1}+[1-h({\lambda }_{k-1}+{\theta }_{k-1})]R_{k-1}\ge 0,\hfill \\ \hfill B_k& =h{\mu }_{k-1}{V}_{k-1}+[1-h({\nu }_{k-1}+{\theta }_{k-1})]B_{k-1}\ge 0.\hfill \end{array}$$

Hence, subtracting I-equation (the third) from A-equation (the second) in (19), we obtain

$$A_k-I_k=h\frac{{\beta }_{k-1}}{N_{k-1}}[S_{k-1}+V_{k-1}]I_k+[1-h({\omega }_{k-1}+{\theta }_{k-1})](A_{k-1}-I_{k-1})\ge 0.$$

On the other hand, $I_{k-1}\le A_{k-1}<N_{k-1}$ and, with the help of (21), we calculate

$$\begin{array}{cc}\hfill I_k& =h{\omega }_{k-1}{A}_{k-1}+[1-h({\gamma }_{k-1}+{\omega }_{k-1}+{\tau }_{k-1}+{\theta }_{k-1})]I_{k-1}\hfill \\ \hfill & <[1-h({\gamma }_{k-1}+{\tau }_{k-1}+{\theta }_{k-1})]N_{k-1}\hfill \\ \hfill & \le N_{k-1}.\hfill \end{array}$$

(22)

Using S-equation (the first) and V-equation (the fifth) in (19) with $t_k-t_{k-1}=h$, we find

$$\begin{array}{cc}\hfill S_k=& h{\Lambda }_{k-1}{N}_{k-1}\hfill \\ & +\left[1-h\left({\alpha }_{k-1}+{\theta }_{k-1}+{\beta }_{k-1}\frac{I_k}{N_{k-1}}\right)\right]S_{k-1}+h{\lambda }_{k-1}{R}_{k-1}+h{\nu }_{k-1}{B}_{k-1},\hfill \\ \hfill V_k=& h{\alpha }_{k-1}{S}_{k-1}+\left[1-h\left({\mu }_{k-1}+{\theta }_{k-1}+{\beta }_{k-1}\frac{I_k}{N_{k-1}}\right)\right]V_{k-1}.\hfill \end{array}$$

Now, using Equations (21) and (22) and since $N_{k-1}-S_{k-1}>0$, we obtain

$$\begin{array}{cc}\hfill S_k\ge & h{\Lambda }_{k-1}(N_{k-1}-S_{k-1})\hfill \\ & +[1+h({\Lambda }_{k-1}-{\alpha }_{k-1}-{\theta }_{k-1}-{\beta }_{k-1})]S_{k-1}+h{\lambda }_{k-1}{R}_{k-1}+h{\nu }_{k-1}{B}_{k-1}>0,\hfill \\ \hfill V_k\ge & h{\alpha }_{k-1}{S}_{k-1}+[1-h({\mu }_{k-1}+{\beta }_{k-1}+{\theta }_{k-1})]V_{k-1}\ge 0.\hfill \end{array}$$

Hence, $N_k>S_k>0.$ The proof is complete. □

Remark 3.

Theorem 3 gives a sufficient condition for non-negativity of the solution to the discrete problem (19). We see that the solution of discrete problem (19) has the same nonnegativity properties as the solution of the continuous problem (17) if the initial data and parameters satisfy analogical conditions, and the step size is small enough (see (21)).

4. Parameter Identification

In a real situation, like the COVID-19 pandemic, the step size in the discrete model (19) is fixed, and we have data for some components of the solution and values of certain parameters. This leads to an inverse problem of a special kind.

In order to formulate a realistic “inverse” problem, that arises from practice, we divide the parameters $p_k$ into two groups

$$p_k=({\tilde{p}}_k,{\widehat{p}}_k),k=1,2,\cdots ,K-1,$$

where ${\tilde{p}}_k:=({\Lambda }_k,{\theta }_k,{\omega }_k,{\mu }_k,{\nu }_k,{\lambda }_k)$ are parameter’s values that it is natural to assume to be known and ${\widehat{p}}_k:=({\alpha }_k,{\beta }_k,{\gamma }_k,{\tau }_k)$ are values that we have to find.

Now, we introduce three groups that do not appear explicitly in SEIRS-BV model, but the official data sets usually contain information for the numbers of individuals in these compartments:

(1)

$Rtotal\left(t\right)$—The cumulative number of the individuals recovered from the disease to the time t. Unlike $R\left(t\right)$, individuals who have already lost disease-acquired immunity are counted in $Rtotal\left(t\right)$.

(2)

$Dtotal\left(t\right)$—The cumulative number of COVID-19 deaths;

(3)

$Vtotal\left(t\right)$—The cumulative number of the fully vaccinated individuals.

It is clear that $Rtotal\left(t\right),Dtotal\left(t\right),andVtotal\left(t\right)$ are not decreasing functions. Actually, in the classical SIR/SEIR models, because of the internal immunity and the absence of vital dynamics, the groups $R\left(t\right)$ and $Rtotal\left(t\right)$ coincide. As in SIR-type models, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{dRtotal}{dt}}& =\gamma \left(t\right)I\left(t\right).\hfill \end{array}$$

(23)

In the same manner (see [25] and references therein),

$$\begin{array}{cc}\hfill {\displaystyle \frac{dDtotal}{dt}}& =\tau \left(t\right)I\left(t\right).\hfill \end{array}$$

(24)

Let us introduce the notation for the available measurements

$$m_k:=(A_k,Rtota{l}_k,Dtota{l}_k,Vtota{l}_k)=(A\left(t_k\right),Rtotal\left(t_k\right),Dtotal\left(t_k\right),Vtotal\left(t_k\right)),$$

$k=1,2,\cdots ,K.$

Similarly, we denote the values of the unknown functions:

$$g_k:=(S_k,I_k,R_k,V_k,B_k).$$

Now, we are ready to formulate the following appropriate “inverse” problem.

Problem IDP: Using the given data $g_1$, ${\left\{{\tilde{p}}_k\right\}}_{k=1}^{K-1}$, ${\left\{m_k\right\}}_{k=1}^K$, find the values

$\left({\left\{g_k\right\}}_{k=2}^K,{\left\{{\widehat{p}}_{k-1}\right\}}_{k=2}^K\right)$, such that the relations (19) hold.

Since the reported data for COVID-19 is on a daily basis, we propose the following algorithm to solve the IDP problem in the case $h=1day$.

Algorithm for solving IDP ( with $\mathbf{h}=\mathbf{1}$):

• Replacing derivatives in (23) and (24) by finite difference with step $h=1$, it leads to the following relations and notations

  $$\begin{array}{cc}\hfill I_{\gamma ,k-1}& :={\gamma }_{k-1}{I}_{k-1}=Rtota{l}_k-Rtota{l}_{k-1},\hfill \\ \hfill I_{\tau ,k-1}& :={\tau }_{k-1}{I}_{k-1}=Dtota{l}_k-Dtota{l}_{k-1},\hfill \end{array}$$

  (25)
  $k=2,3,\cdots ,K$. Since the values ${\left\{Dtotal{}_k\right\}}_{k=1}^K$ and ${\left\{Rtotal{}_k\right\}}_{k=1}^K$ are given and nondecreasing in respect to k, we find via (25) the nonnegative values ${\left\{I_{\gamma ,k-1}\right\}}_{k=2}^K$, ${\left\{I_{\tau ,k-1}\right\}}_{k=2}^K$.
• Since $N_1=S_1+A_1+R_1+V_1+B_1$ is given, the relations (20), (24) and (25) imply

  $$N_k=(1+{\Lambda }_{k-1}-{\theta }_{k-1})N_{k-1}-I_{\tau ,k-1},k=2,3,\cdots ,K.$$

  (26)
  Thus the values of population size can be calculated.
• The values of the vaccination parameter are

  $${\alpha }_{k-1}=\sigma \frac{Vtota{l}_k-Vtota{l}_{k-1}}{N_{k-1}},k=2,3,\cdots ,K,$$

  (27)

  where $\sigma$ is the vaccine effectiveness. Since the values ${\left\{Vtotal{}_k\right\}}_{k=1}^K$ are given, we find via (27) the nonnegative values ${\left\{{\alpha }_{k-1}\right\}}_{k=2}^K$.
• Obviously, $N_{k-1}\ne 0$ and we introduce the notations

  $$I_{\beta ,k-1}:={\beta }_{k-1}\frac{I_k}{N_{k-1}},k=2,3,\cdots ,K.$$

  (28)
  Now, going step by step from $k-1$ to k, using (25) and (28), we rewrite the relations (19) in the form

  $$\left|\begin{array}{c}{S}_k=[1-{\alpha }_{k-1}-{\theta }_{k-1}-I_{\beta ,k-1}]S_{k-1}+{\Lambda }_{k-1}{N}_{k-1}+{\lambda }_{k-1}{R}_{k-1}+{\nu }_{k-1}{B}_{k-1},\hfill \\ \\ A_k=(1-{\theta }_{k-1})A_{k-1}+(S_{k-1}+V_{k-1})I_{\beta ,k-1}-I_{\gamma ,k-1}-I_{\tau ,k-1},\hfill \\ \\ I_k=[1-{\omega }_{k-1}-{\theta }_{k-1}]I_{k-1}+{\omega }_{k-1}{A}_{k-1}-I_{\gamma ,k-1}-I_{\tau ,k-1}\hfill \\ \\ R_k=[1-{\lambda }_{k-1}-{\theta }_{k-1}]R_{k-1}+I_{\gamma ,k-1},\hfill \\ \\ V_k=[1-{\mu }_{k-1}-{\theta }_{k-1}-I_{\beta ,k-1}]V_{k-1}+{\alpha }_{k-1}{S}_{k-1},\hfill \\ \\ B_k=[1-{\nu }_{k-1}-{\theta }_{k-1}]B_{k-1}+{\mu }_{k-1}{V}_{k-1},\hfill \end{array}\right.$$

  (29)

  where $k=2,3,\dots ,K$.
  Starting with the nonnegative initial data $A_1$ and $g_1=(S_1,I_1,R_1,V_1,B_1),$ where $S_1>0$ and $A_1\ge I_1>0$, we calculate

  4.1.

  From the second equation in (29), we obtain

  $$I_{\beta ,1}=\frac{1}{S_1+V_1}[A_2+({\theta }_1-1)A_1+I_{\gamma ,1}+I_{\tau ,1}].$$

  We note that $I_{\beta ,1}\ge 0$, because $A_2-A_1+I_{\gamma ,1}+I_{\tau ,1}\ge 0$ are the new cases for day $t_2$.

  4.2.

  Now, using (29), we calculate consistently the values ${\left\{S_k\right\}}_{k=2}^K$, ${\left\{I_k\right\}}_{k=2}^K$, ${\left\{R_k\right\}}_{k=2}^K$, ${\left\{V_k\right\}}_{k=2}^K$, ${\left\{B_k\right\}}_{k=2}^K$ and, if $S_{k-1}+V_{k-1}\ne 0$

  $$I_{\beta ,k-1}=\frac{1}{S_{k-1}+V_{k-1}}[A_k+({\theta }_{k-1}-1)A_{k-1}+I_{\gamma ,k-1}+I_{\tau ,k-1}],k=3,4,\cdots ,K.$$

  (30)

  4.3.

  Finally, if $I_{k-1}\ne 0$ and $I_k\ne 0$, we calculate

  $${\beta }_{k-1}=N_{k-1}\frac{I_{\beta ,k-1}}{I_k},{\gamma }_{k-1}=\frac{I_{\gamma ,k-1}}{I_{k-1}},{\tau }_{k-1}=\frac{I_{\tau ,k-1}}{I_{k-1}},k=2,3,\cdots ,K.$$

  (31)

  4.4.

  If one or more of the values $S_{k-1}+V_{k-1}$, $I_{k-1}$ is equal to zero for some k, and the algorithm must stop at the first such k. Otherwise, the algorithm continues and all values ${\left\{g_k\right\}}_{k=1}^K$ and ${\left\{{\widehat{p}}_{k-1}\right\}}_{k=2}^K$ can be uniquely determined. Actually, since we are looking for a biologically reasonable solution of the problem $IDP$, i.e., for nonnegative solution of this problem, we should also stop the algorithm if some of the calculated values in the step 4.2 are negative.

5. Identification of Parameters in the Discrete Problem

We conduct a set of numerical experiments to solve $IDP$ problem with Bulgarian COVID-19 data and to study the original differential problem (1). The period under consideration is 8 March 2020–10 April 2022. Within the above-mentioned period, we extract several parts that correspond to domination of different variants of the coronavirus. In each of them, we fix the values of the known parameters ${\left\{{\tilde{p}}_k\right\}}_{k=1}^{K-1}$ using some official data respectively: [33,34,35,36,37,38]. In order to calculate appropriate values for parameters ${\tilde{p}}_k:=({\Lambda }_k,{\theta }_k,{\omega }_k,{\mu }_k,{\nu }_k,{\lambda }_k)$, we suppose that

• ${\Lambda }_k=\Lambda$ is the average birth rate for 2015–2020;
• ${\theta }_k=\theta$ is the average natural mortality rate for 2015–2020;
• ${\omega }_k=1/T_e$, where $T_e$ is the incubation (latency) period for the dominant variant of COVID-19;
• ${\mu }_k=1/T_a$, where $T_a$ is the average time taken for antibodies to develop;
• ${\nu }_k=1/T_b$, where $T_b$ is the duration of the immune responses in individuals with vaccination-acquired immunity;
• ${\lambda }_k=1/T_r$, where $T_r$ is the duration of the immune responses in recovered individuals.

For calculation of the birth and natural mortality rates, we use available data in the Information System INFOSTAT of the National Statistical Institute of the Republic of Bulgaria [35].

Four vaccines are in use in Bulgaria—Comirnaty (developed by BioNTech and Pfizer) known as Pfizer, Vaxzevria (previously COVID-19 Vaccine AstraZeneca), Spikevax (previously COVID-19 Vaccine Moderna) and Janssen. Taking into account the numbers of the fully vaccinated people with these vaccines and the product information in [38], we calculate the average values of vaccines’ parameters. This way, we can assume that, during the COVID-19 mass vaccination campaign in Bulgaria, one vaccine is used with average values of parameters specified in

Table 2. The given parameter’s values.

Parameter	Description	Values
$\Lambda$	birth rate	$2.4095\times {10}^{-5}$
$\theta$	natural mortality rate	$4.1904\times {10}^{-5}$
$T_e$	latency period	7 days 1, 6 days 2, 5 days 3, 4 days 4
$T_a$	time taken for antibodies to develop	14 days
$T_b$	duration of vaccine-based immunity	180 days
$T_r$	duration of disease-based immunity	180 days
$\sigma$	vaccine effectiveness	${0.85}^2,{0.70}^3,{0.45}^4$

for ¹ Wuhan variant, ² Alpha variant, ³ Delta variant, ⁴ Omicron variant.

.

The initial data for COVID-19 pandemic in Bulgaria are

$$A_1=4,g_1=(S_1,I_1,R_1,V_1,B_1)=(6941259,4,0,0,0).$$

(32)

We apply the algorithm described in Section 4 with step size $h=1$, officially reported measurements ${\left\{m_k\right\}}_{k=1}^K={\left\{(A_k,Rtota{l}_k,Dtota{l}_k,Vtota{l}_k)\right\}}_{k=1}^K$, initial data (32), and the parameters’ values ${\left\{{\tilde{p}}_k\right\}}_{k=1}^{K-1}$ (Table 2). The algorithm provides a unique non-negative solution

$${\left\{g_k\right\}}_{k=1}^K={\left\{(S_k,I_k,R_k,V_k,B_k)\right\}}_{k=1}^K,{\left\{{\widehat{p}}_{k-1}\right\}}_{k=2}^K={\left\{({\alpha }_{k-1},{\beta }_{k-1},{\gamma }_{k-1},{\tau }_{k-1})\right\}}_{k=2}^K$$

of problem $IDP$ and $S_k>0,A_k\ge I_k>0$ for $k=1,2,\cdots ,K$. The obtained weekly average values of infection and recovery rates are given in Figure 2.

Remark 4.

The Wuhan is usually used as a name of the first variant of the virus. To avoid possible misunderstanding, we note that, in the following Figure 2, Figure 3 and Figure 4, we point out the names of different variants (Wuhan, Alpha, Delta, Omicron) of SARS-CoV-2 and the corresponding periods in which they dominated in Bulgaria.

We observe that the infection rate ${\beta }_{week}$ decreases at the time of partial lockdowns (27 October 2020, 22 March 2021) or when strong rules and restrictions (21 October 2021) were imposed. The above-mentioned infection rate has been decreasing since green certificates were introduced (21 October 2021) as a preventive measure. On the other hand, the infection rate increased rapidly with the beginning of the first (the autumn of 2020) and second (the spring of 2021) waves and with the invasion of Delta (22 July 2021) or Omicron (2 January 2022) variants of SARS-CoV-2 in Bulgaria. We observe the following relations for the infection rate values for the first two weeks of Delta ($\delta$) or Omicron (o) waves ${\beta }_{o,wee{k}_1}\approx 2{\beta }_{\delta ,wee{k}_1}$ and ${\beta }_{o,wee{k}_2}\approx 3{\beta }_{\delta ,wee{k}_2}$. The variation in infection rate for the Omicron variant is less than for the other virus variants. The reason for this is that the measures of social distancing were weak during the Omicron wave. This behavior of the Omicron infection rate causes a big wave which reaches its peak faster than the Delta wave. At the same time, the peaks of the recovery rate ${\gamma }_{week}$ follow the peaks of the infection rate with some delay.

The obtained weekly average values of the vaccination parameter ${\alpha }_{week}$ and mortality rate ${\tau }_{week}$ are given in Figure 3.

We observe that the vaccination parameter ${\alpha }_{week}$ increases with the introduction of green certificates in Bulgaria, but it is still the lowest rate in the European Union and one of the lowest rates globally. At the same time, the COVID-19 mortality rate in Bulgaria is one of the highest levels in the world. It is natural that the peaks of the mortality rate ${\tau }_{week}$ follow the peaks of the infection rate with some delay.

6. Numerical Solution of the Differential Problem

Using Bulgarian COVID-19 data, in Section 5, we already found the daily values ${\alpha }_k,{\beta }_k,{\gamma }_k,{\tau }_k$, $k=1,2,\cdots ,K-1$ of the parameters in the discrete problem (19). A natural question that arises is: How close are they to the daily values of the parameters $\alpha \left(t\right),\beta \left(t\right),\gamma \left(t\right),\tau \left(t\right)$ in the differential problem (1)?

In order to answer this question, we solve numerically a recurrent sequence of initial problems for differential equation systems. Each of these problems corresponds to one of the days $t_1,t_2,\dots ,t_K$ of the considered period of the pandemic. We do this using the selected values of parameters ${\left\{{\Lambda }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\theta }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\omega }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\lambda }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\nu }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\mu }_k\right\}}_{k=1}^{K-1}$ and the calculated values for parameters ${\left\{{\alpha }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\beta }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\gamma }_k\right\}}_{k=1}^{K-1}$, ${\left\{{\tau }_k\right\}}_{k=1}^{K-1}$. More precisely, using the built-in function ode45 in Matlab (see [39]), which implements a version of Runge–Kutta 4th/5th-order method, for $t\in [t_{k-1},t_k]$, we solve the Cauchy problem $SEIRS-V{B}_k$, $k=2,3,\dots ,K$:

$$\left|\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{S}}_k}{dt}=}& {\Lambda }_{k-1}{\tilde{N}}_k\left(t\right)-({\alpha }_{k-1}+{\theta }_{k-1}){\tilde{S}}_k\left(t\right)-\frac{{\beta }_{k-1}}{{\tilde{N}}_k\left(t\right)}{\tilde{S}}_k\left(t\right){\tilde{I}}_k\left(t\right)\hfill \\ \hfill & +{\lambda }_{k-1}{\tilde{R}}_k\left(t\right)+{\nu }_{k-1}{\tilde{B}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{E}}_k}{dt}=}& \frac{{\beta }_{k-1}}{{\tilde{N}}_k\left(t\right)}{\tilde{S}}_k\left(t\right){\tilde{I}}_k\left(t\right)+\frac{{\beta }_{k-1}}{{\tilde{N}}_k\left(t\right)}{\tilde{V}}_k\left(t\right){\tilde{I}}_k\left(t\right)-({\omega }_{k-1}+{\theta }_{k-1}){\tilde{E}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{I}}_k}{dt}=}& {\omega }_{k-1}{\tilde{E}}_k\left(t\right)-[{\gamma }_{k-1}+{\tau }_{k-1}+{\theta }_{k-1}]{\tilde{I}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{R}}_k}{dt}=}& {\gamma }_{k-1}{\tilde{I}}_k\left(t\right)-({\lambda }_{k-1}+{\theta }_{k-1}){\tilde{R}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{V}}_k}{dt}=}& {\alpha }_{k-1}{\tilde{S}}_k-\frac{{\beta }_{k-1}}{{\tilde{N}}_k\left(t\right)}{\tilde{V}}_k\left(t\right){\tilde{I}}_k\left(t\right)-({\mu }_{k-1}+{\theta }_{k-1}){\tilde{V}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{B}}_k}{dt}=}& {\mu }_{k-1}{\tilde{V}}_k\left(t\right)-({\nu }_{k-1}+{\theta }_{k-1}){\tilde{B}}_k\left(t\right),\hfill \\ \\ \hfill {\displaystyle \frac{d{\tilde{N}}_k}{dt}=}& ({\Lambda }_{k-1}-{\theta }_{k-1}){\tilde{N}}_k\left(t\right)-{\tau }_{k-1}{\tilde{I}}_k\left(t\right),\hfill \\ \\ \hfill {\tilde{S}}_k\left(t_{k-1}\right)=& {\tilde{S}}_{k-1}\left(t_{k-1}\right),{\tilde{E}}_k\left(t_{k-1}\right)={\tilde{E}}_{k-1}\left(t_{k-1}\right),{\tilde{I}}_k\left(t_{k-1}\right)={\tilde{I}}_{k-1}\left(t_{k-1}\right),\hfill \\ \hfill {\tilde{R}}_k\left(t_{k-1}\right)=& {\tilde{R}}_{k-1}\left(t_{k-1}\right),{\tilde{V}}_k\left(t_{k-1}\right)={\tilde{V}}_{k-1}\left(t_{k-1}\right),{\tilde{B}}_k\left(t_{k-1}\right)={\tilde{B}}_{k-1}\left(t_{k-1}\right),\hfill \\ \hfill {\tilde{N}}_k\left(t_{k-1}\right)=& {\tilde{N}}_{k-1}\left(t_{k-1}\right),\hfill \end{array}\hfill \end{array}\right.$$

(33)

where ${\tilde{S}}_1\left(t_1\right)=S_1,{\tilde{E}}_1\left(t_1\right)=E_1=0,{\tilde{I}}_1\left(t_1\right)=I_1=A_1,$ ${\tilde{R}}_1\left(t_1\right)=R_1,{\tilde{V}}_1\left(t_1\right)=V_1,{\tilde{B}}_1\left(t_1\right)=B_1,{\tilde{N}}_1\left(t_1\right)=N_1=S_1+E_1+I_1+R_1+V_1+B_1.$

The values of functions ${\tilde{S}}_{k-1}\left(t\right),$${\tilde{E}}_{k-1}\left(t\right),$ ${\tilde{I}}_{k-1}\left(t\right),$${\tilde{R}}_{k-1}\left(t\right),$ ${\tilde{V}}_{k-1}\left(t\right)$, ${\tilde{B}}_{k-1}\left(t\right)$ and ${\tilde{N}}_{k-1}\left(t\right)$ (solution of the problem $SEIRS-V{B}_{k-1}$) at the end of the day $t_{k-1}$ are used as initial data in the next problem $SEIRS-V{B}_k,$ which describe the quantitative changes over the day $t_k.$ The obtained piecewise linear curve is shown in Figure 4.

Since the numbers of active cases $A_{pandemic}=(A_1,A_2,\cdots ,A_K)$ are known, we will compare them with the values ${\tilde{A}}_{pandemic}=({\tilde{E}}_1\left(t_1\right)+{\tilde{I}}_1\left(t_1\right),{\tilde{E}}_2\left(t_2\right)+{\tilde{I}}_2\left(t_2\right),\cdots ,{\tilde{E}}_K\left(t_K\right)+{\tilde{I}}_K\left(t_K\right))$, obtained by the procedure, described above. To do this, we use two different norms—the ${\ell }_2$ one and the sup-norm ${\ell }_{\infty }$, i.e.,

$$Error({\ell }_2,A):=\frac{\parallel {\tilde{A}}_{pandemic}-A_{pandemic}{\parallel }_2}{\parallel A_{pandemic}{\parallel }_2}=0.044,$$

(34)

$$Error({\ell }_{\infty },A):=\frac{\parallel {\tilde{A}}_{pandemic}-A_{pandemic}{\parallel }_{\infty }}{\parallel A_{pandemic}{\parallel }_{\infty }}=0.035.$$

(35)

These results show that the suggested method for identification of the parameters in the discrete problem leads to the parameters’ values which are very close to the values of the time-dependent coefficients in the differential problem. We can use other discretizations of the SEIRS-VB model, based for example on an explicit or implicit Euler method, but better results give the suggested semi-implicit discretized model (19).

7. Short-Term Forecasting

In this section, we present a methodology for short-term prediction of the numbers of active cases (A), new daily cases, COVID-19 deaths ($Dtotal$) and cumulative number or recovered individuals ($Rtotal$). We assume that the parameters’ values ${\tilde{p}}_k:=({\Lambda }_k,{\theta }_k,{\omega }_k,{\mu }_k,{\nu }_k,{\lambda }_k)$ are known. The experiments are based on the parameter identification method developed in Section 4 and the results obtained in Section 5. One possible way is to predict the unknown parameters’ values ${\widehat{p}}_k:=({\alpha }_k,{\beta }_k,{\gamma }_k,{\tau }_k)$ for the short- term forecast period using their values (calculated in Section 5) for several days before this period. It is worth mentioning that, in order to solve the problem $IDP$ (see Algorithm for solving IDP in Section 4), we have to calculate the values

$$u_{k-1}:=\left(I_{\beta ,k-1},I_{\gamma ,k-1},I_{\tau ,k-1},{\alpha }_{k-1}\right),$$

first and then, via (31), we obtain ${\beta }_{k-1},{\gamma }_{k-1},{\tau }_{k-1}.$ For that purpose, we predict the values of $u_{k-1}$ in order to make the forecast for numbers of individuals in the compartments of SEIRS-VB model more accurate. To fix the notations, let the first day of the forecast period be denoted by $t_n$. Then, the last day of this forecast period will be denoted $t_{n+13}$.

We denote by

$${\tilde{u}}_{k-1}:=\left({\check{I}}_{\beta ,k-1},{\check{I}}_{\gamma ,k-1},{\check{I}}_{\tau ,k-1},{\check{\alpha }}_{k-1}\right)$$

the predicted values of the parameters, using the proposed methodology, in order to distinguish them from the corresponding official ground-true data $u_{k-1}$, reported by the Bulgarian officials.

The last known values are $A_{n-1},Rtota{l}_{n-1},Dtota{l}_{n-1}$ and the last calculated values are $S_{n-1},I_{n-1},R_{n-1},V_{n-1},B_{n-1}$, which are obtained using the last known parameter’s value $u_{n-2}$. We will predict the values ${\check{S}}_k,{\check{A}}_k,{\check{I}}_k,{\check{R}}_k,{\check{V}}_k,{\check{B}}_k$, $\check{R}tota{l}_k,\check{D}tota{l}_k$ and ${\check{u}}_{k-1}$ for $k=n,n+1,\cdots ,n+13$.

We perform two sets of numerical experiments, related to predicting a 14-day-time-frame in the future, based on the available official data for Bulgaria up to the first day of the time-frame. These time-frames are: 18–31 October 2021 and 7–20 February 2022, respectively. During the first time-frame, the Delta variant of the virus was dominant, while, during the second time-frame, the Omicron variant dominated.

For each time-frame, we make the prediction in three steps:

Prediction step 1. It is easy to discover cases/weekdays dependency in the official data. For example, the officially reported cases for Saturday and Sunday (reported on Sunday and Monday, respectively) are much fewer than the cases for Monday (reported on Tuesday). The same is true if we compare a holiday with a working day of the week. For this reason, we assume dependency of the parameters’ values on the weekdays. To verify this assumption, we examine the ratios between any two consecutive daily parameter values over different 7-day periods before the forecast period. For example, we study the behavior of the ratios

$$\begin{array}{c}\hfill {\delta }_{\beta ,k}:=({\delta }_{\beta ,k}^1,{\delta }_{\beta ,k}^2,\cdots ,{\delta }_{\beta ,k}^7)=(I_{\beta ,k}/I_{\beta ,k-1},I_{\beta ,k+1}/I_{\beta ,k},\cdots ,I_{\beta ,k+6}/I_{\beta ,k+5})\end{array}$$

(36)

and ${\delta }_{\gamma ,k},{\delta }_{\tau ,k}$, ${\delta }_{\alpha ,k},$ defined in a similar way for $\gamma ,\tau ,\alpha$, respectively. Let us note that all these ratios are known for $k\le n-8$. After analyzing them, we forecast in an appropriate way the values ${\check{\delta }}_{\beta ,n-1},{\check{\delta }}_{\beta ,n+6}$, ${\check{\delta }}_{\gamma ,n-1},{\check{\delta }}_{\gamma ,n+6}$, ${\check{\delta }}_{\tau ,n-1},{\check{\delta }}_{\tau ,n+6}$, ${\check{\delta }}_{\alpha ,n-1},{\check{\delta }}_{\alpha ,n+6}$.

Prediction step 2. The first predicted day is $t_n$ and the values of parameters $I_{\beta ,n-2}$$I_{\gamma ,n-2},I_{\tau ,n-2},{\alpha }_{n-2}$ for the day $t_{n-1}$ are known. Then, our forecast of parameters’ values for the day $t_n$ is

$$\begin{array}{cc}\hfill & {\check{I}}_{\beta ,n-1}={\check{\delta }}_{\beta ,n-1}^1{I}_{\beta ,n-2},{\check{I}}_{\gamma ,n-1}={\check{\delta }}_{\gamma ,n-1}^1{I}_{\gamma ,n-2},\hfill \\ \hfill & {\check{I}}_{\tau ,n-1}={\check{\delta }}_{\tau ,n-1}^1{I}_{\tau ,n-2},{\check{\alpha }}_{n-1}={\check{\delta }}_{\alpha ,n-1}^1{\alpha }_{n-2}\hfill \end{array}$$

(37)

and for the next six days are as follows:

$$\begin{array}{cc}\hfill & {\check{I}}_{\beta ,n+j-2}={\check{\delta }}_{\beta ,n-1}^j{\check{I}}_{\beta ,n+j-3},{\check{I}}_{\gamma ,n+j-2}={\check{\delta }}_{\gamma ,n-1}^j{\check{I}}_{\gamma ,n+j-3},\hfill \\ \hfill & {\check{I}}_{\tau ,n+j-2}={\check{\delta }}_{\tau ,n-1}^j{\check{I}}_{\tau ,n+j-3},{\alpha }_{n+j-2}={\check{\delta }}_{\alpha ,n-1}^j{\check{\alpha }}_{n+j-3}\hfill \end{array}$$

(38)

for $j=2,3,\cdots ,7$. Similarly, for the next 7 days, we define

$$\begin{array}{cc}\hfill & {\check{I}}_{\beta ,n+6}={\check{\delta }}_{\beta ,n+6}^j{\check{I}}_{\beta ,n+j+4},{\check{I}}_{\gamma ,n+j+5}={\check{\delta }}_{\gamma ,n+6}^j{\check{I}}_{\gamma ,n+j+4},\hfill \\ \hfill & {\check{I}}_{\tau ,n+j+5}={\check{\delta }}_{\tau ,n+6}^j{\check{I}}_{\tau ,n+j+4},{\alpha }_{n+j+5}={\check{\delta }}_{\alpha ,n+6}^j{\check{\alpha }}_{n+j+4}\hfill \end{array}$$

(39)

for $j=1,2,\cdots ,7$.

In such a way, we predict all parameters’ values which are needed for calculating the number of individuals in each of the considered compartments.

Prediction step 3. Since the values $S_{n-1},A_{n-1},I_{n-1},R_{n-1},V_{n-1},B_{n-1}$ and $Rtota{l}_{n-1}$, $Dtota{l}_{n-1}$ for the last day before the considered forecast period are known, we can predict the number of individuals in each of these compartments during the period $t_n,\cdots ,t_{n+13}$. Using (25) and the predicted daily values of parameters, we calculate the values ${\left\{\check{R}toata{l}_k\right\}}_{k=n}^{n+13}$ and ${\left\{\check{D}toata{l}_k\right\}}_{k=n}^{n+13}$ of the recovered individuals and COVID-19 deaths, respectively. Furthermore, using again the predicted daily values of the parameters and (29), we calculate the values ${\left\{{\check{A}}_k\right\}}_{k=n}^{n+13}$ of active cases and the new daily cases ${\{{\check{A}}_k-{\check{A}}_{k-1}+{\check{I}}_{\gamma ,k-1}+{\check{I}}_{\tau ,k-1}\}}_{k=n+1}^{n+13}$. More precisely, for each day of the considered time-frame, using the morning values $({\check{S}}_{k-1},{\check{A}}_{k-1},{\check{I}}_{k-1},{\check{R}}_{k-1},{\check{V}}_{k-1},{\check{B}}_{k-1})$, $(\check{D}tota{l}_{k-1},\check{R}tota{l}_{k-1})$, the parameter values ${\tilde{p}}_{k-1}$ and the predicted values ${\check{u}}_{k-1}$, we calculate $({\check{S}}_k,{\check{A}}_k,{\check{I}}_k,{\check{R}}_k,{\check{V}}_k,{\check{B}}_k)$ (using (29)) and $(\check{D}tota{l}_k,\check{R}tota{l}_k)$ (using (25)).

7.1. The First Time-Frame 18–31 October 2021

This time-frame is related to the peak of new daily cases of the wave caused by Delta variant of the virus. The reported cases on October 19th are considered as the values for October 18 th and so on. We consider the known ratios ${\delta }_{\beta ,k},{\delta }_{\gamma ,k},{\delta }_{\tau ,k}$ and ${\delta }_{\alpha ,k}$ for the last seven weeks before the first predicted day—18th October 2021, Monday. The values of the ratio ${\delta }_{\beta ,k}$ are given in

Table 3. Documentation of the ratio ${\delta }_{\beta ,k}$.

	Mon/Sun	Tue/Mon	Wed/Tue	Thur/Wed	Fri/Thur	Sat/Fri	Sun/Sat
k	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^1$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^2$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^3$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^4$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^5$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^6$	${\mathit{\delta }}_{\mathit{\beta },\mathit{k}}^7$
$n-8$	5.573	0.999	0.959	0.957	1.081	0.543	0.536
$n-15$	5.452	1.021	0.887	0.905	1.079	0.502	0.499
$n-22$	5.098	1.014	0.832	0.957	1.020	0.587	0.418
$n-29$	5.138	0.929	0.388	2.156	1.037	0.545	0.478
$n-36$	4.559	1.015	0.918	0.926	1.097	0.503	0.482
$n-43$	1.386	2.957	0.975	0.833	0.998	0.559	0.481
$n-50$	4.624	0.982	0.850	0.896	1.048	0.521	0.572

.

As we expected, almost always ${\delta }_{\beta ,k}^1>{\delta }_{\beta ,k}^j,j=2,3,\cdots ,7$. Of course, there is an exception: ${\delta }_{\beta ,n-43}^1<{\delta }_{\beta ,n-43}^2$ because Monday, 6th September 2021 (Unification Day) is an official holiday in Bulgaria, and the reported new infected cases are much fewer than usual. An analogous situation is with ${\delta }_{\beta ,n-29}^3$ and ${\delta }_{\beta ,n-29}^4$ because 22nd September 2021 (The Independence Day in Bulgaria) is also a holiday. Such values should not be used in prediction of ${\check{\delta }}_{\beta ,n-1}$ and ${\check{\delta }}_{\beta ,n+6}$.

On the other hand, the values in Table 3 show that, in each column, ${\delta }_{\beta ,k}^j\approx cons{t}_j$ during almost all considered weeks. That is why we set in prediction step 1

$$\begin{array}{cc}\hfill {\check{\delta }}_{\beta ,n-1}^j& :=\frac{1}{3}\left({\delta }_{\beta ,n-8}^j+{\delta }_{\beta ,n-15}^j+{\delta }_{\beta ,n-22}^j\right),\hfill \\ \hfill {\check{\delta }}_{\beta ,n+6}^j& :=\frac{1}{3}\left({\check{\delta }}_{\beta ,n-1}^j+{\delta }_{\beta ,n-8}^j+{\delta }_{\beta ,n-15}^j\right)\hfill \end{array}$$

(40)

for $j=1,2,\cdots ,7.$

Now, using (37)–(39) in prediction step 2, we derive the desired values ${\left\{{\check{I}}_{\beta ,k-1}\right\}}_{k=n}^{n+13}$. The behavior of the other ratios in the considered seven weeks is similar, and, in the same manner (as (40)), we predict other daily parameters’ values. Finally, we can make the prediction step 3.

The experimental comparison between the official data and our prediction for the first time-frame is illustrated in Figure 5, Figure 6, Figure 7 and Figure 8. The official data are in the blue color and the predicted values are in the red color.

We observe a very good agreement between the predicted and the confirmed data during the first and the second week of the time-frame. However, the forecast for the first week is more accurate than for the second. It is interesting that, with this experiment, we predict the exact day of the peak of new daily cases and even the numbers of the predicted and reported cases on this day are very close. The number of recovered individuals and COVID-19 deaths are also fitted very well.

In order to analyze the accuracy of the proposed prediction methodology, the relative weekly-based errors (see (34) and (35)) for the compartments for both the first and the second 14-days-time-frames are documented in

Table 4. Relative errors of the prediction for week 1 (18–24 October 2021) and week 2 (25–31 October 2021).

Compartment	Error (${\mathit{l}}_2$, Week${}_1$)	Error (${\mathit{l}}_2$, Week${}_2$)	Error (${\mathit{l}}_{\mathit{\infty }}$, Week${}_1$)	Error (${\mathit{l}}_{\mathit{\infty }}$, Week${}_2$)
Active cases	0.018	0.034	0.029	0.042
New daily cases	0.090	0.138	0.132	0.168
Deaths	0.005	0.013	0.007	0.016
Recovered	0.001	0.001	0.002	0.002

.

7.2. The Second Time-Frame 7–20 February 2022

This time-frame is related to the peak of active cases of the wave caused by the Omicron variant of the virus. Actually, in Bulgaria, this is the pandemic’s tallest peak (see Figure 4). Because of the very fast invasion of the Omicron variant, the peak has been reached only four weeks after the first confirmed cases with this variant. We consider again the known ratios ${\delta }_{\beta ,k},{\delta }_{\gamma ,k},{\delta }_{\tau ,k}$ and ${\delta }_{\alpha ,k}$ for the last seven weeks before the first predicted day—7th February 2022. Behavior of ${\delta }_{\beta ,k},{\delta }_{\tau ,k},{\delta }_{\alpha ,k}$ is similar to the Delta wave, wile ${\delta }_{\gamma ,k}$ oscillates (see

Table 5. Documentation of the ratio ${\delta }_{\gamma ,k}$.

	Mon/Sun	Tue/Mon	Wed/Tue	Thur/Wed	Fri/Thur	Sat/Fri	Sun/Sat
k	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^1$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^2$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^3$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^4$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^5$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^6$	${\mathit{\delta }}_{\mathit{\gamma },\mathit{k}}^7$
$n-8$	2.311	0.857	1.051	1.306	1.135	0.557	0.593
$n-15$	4.249	0.459	1.790	0.971	0.656	0.292	2.444
$n-22$	0.447	3.186	0.490	0.442	5.049	0.649	1.101
$n-29$	10.917	1.425	4.729	0.650	0.562	0.408	1.326
$n-36$	1.229	2.100	0.565	2.221	0.401	2.655	0.047
$n-43$	4.265	1.381	1.485	0.485	0.512	0.201	4.522
$n-50$	8.402	0.878	0.734	0.950	0.237	0.681	1.293

).

Since the behavior of ${\delta }_{\gamma ,n-22}^j$ is very different from behavior of ${\delta }_{\gamma ,n-8}^j$ and ${\delta }_{\gamma ,n-15}^j$ in contrast to the first time-frame (see (40)) in prediction step 1, we set

$${\check{\delta }}_{\gamma ,n-1}^j:=\frac{1}{2}\left({\delta }_{\gamma ,n-8}^j+{\delta }_{\gamma ,n-15}^j\right),{\check{\delta }}_{\gamma ,n+6}^j:=\frac{1}{2}\left({\check{\delta }}_{\gamma ,n-1}^j+{\delta }_{\gamma ,n-8}^j\right),j=1,2,\cdots ,7.$$

In a similar way, we calculate ${\check{\delta }}_{\beta ,n-1}^j,{\check{\delta }}_{\tau ,n-1}^j,{\check{\delta }}_{\alpha ,n-1}^j$ and ${\check{\delta }}_{\beta ,n+1}^j,{\check{\delta }}_{\tau ,n+1}^j,{\check{\delta }}_{\alpha ,n+6}^j$ for $j=1,2,\cdots ,7.$

To predict the daily values of parameters (prediction step 2), we use the same methodology as in the prediction of the first time-frame. Now, we are able to make a prediction step 3. The comparison between the official data and our prediction for the second time-frame is illustrated in Figure 9, Figure 10, Figure 11 and Figure 12. The official data are in the blue color and the predicted data are in the red color.

We observe a good agreement between the predicted and the reported data. The peak of active cases and the number of COVID-19 deaths are fitted very well. The predicted active cases for the second week are smaller than the officially documented ones and the daily differences between the two plots are larger than in the first time-frame. At the same time, the predicted number of new cases is larger than the reported ones. The number of recovered individuals is fitted well. This experiment shows that the proposed methodology works well even in such a delicate situation, when the parameters change their behavior quite significantly.

The relative weekly-based errors for the compartments for both the first and the second 14-days-time-frames are documented in

Table 6. Relative errors of the prediction for week 1 (7–13 February 2022) and week 2 (14–20 February 2022).

Compartment	Error (${\mathit{l}}_2$, Week${}_1$)	Error (${\mathit{l}}_2$, Week${}_2$)	Error (${\mathit{l}}_{\mathit{\infty }}$, Week${}_1$)	Error (${\mathit{l}}_{\mathit{\infty }}$, Week${}_2$)
Active cases	0.015	0.048	0.023	0.077
New daily cases	0.107	0.230	0.137	0.210
Deaths	0.001	0.005	0.001	0.007
Recovered	0.001	0.028	0.002	0.041

.

We observe a good agreement between the corresponding error margins of the two norms $l_2$ and $l_{\infty }$ (see (34) and (35), Table 4 and Table 6), which suggest robustness of the proposed model with respect to the norm choice within the $\left\{l_p\right\}$ family. Let us note that, for the active cases, deaths and recovered individuals, we predict the cumulative number, while, for the dally cases, we consider only the new infectious individuals on the corresponding day. Thus, the relative forecast errors are naturally larger for the new daily cases than the corresponding forecast errors for the other compartments.

8. Discussion

Actually, since the considered model takes into account some factors that play an essential role in the viral diseases’ dynamics such as re-susceptibility, duration of the latency period (specific for each of the dominant virus variants) and the vaccines’ effectiveness. After finding the models parameters in Section 5, we are able to calculate the daily number of Susceptible individuals $S_k$, Recovered individuals $R_k$, Vaccinated susceptible individuals $V_k$ and individuals with antibodies $B_k$. It should be noted that there is no official data on the size of these compartments, as public data sets contain information on the cumulative number of recovered and vaccinated persons. That is why the results obtained by our model can be very useful. For example, in Figure 13, we give the daily relative size of the compartment of all susceptible individuals $S_k+V_k$ in Bulgaria, which is closely related to the determination of the herd immunity threshold (see [40]).

In the months following the emergence of SARS-CoV-2, “herd immunity” was frequently cited as the long-term destination of the COVID-19 pandemic. However, since the “Delta” variant appeared (end of August 2021), our idea was to calculate in our mathematical-computer model the level of “herd immunity”, and it seemed to be more than 90%, which is impossible to be achieved. From our group, we announced this fact and discussed it at a Conference of UBM on 5 September 2021. At the same time, at least two publications [41,42] in this area appeared in the USA with discussions on how it is better to react to this situation?

The idea of herd immunity primarily supports high vaccination coverage and the acquisition of natural immunity due to illness.

Although COVID-19 vaccines provide some protection against infection and a mild form of COVID-19, they have failed to stop the transmission of the virus, especially for the highly transmitted delta and omicron variants. An excellent example in this regard was the decision of the United Kingdom Government to allow the opening of society based on a high percentage of vaccination coverage, which minimizes the risk of severe COVID-19 and death. Indeed, significantly increased mortality was not achieved, but the daily incidence reached over 270,000 new cases per day with Omicron. This decision has demonstrated that the purpose of herd immunity, even in a resource-rich environment, is unattainable [43].

An interesting study based on 17 cases of genetically confirmed re-infection with COVID-19 showed that one immunocompromised patient had mild symptoms in the first infection but developed severe symptoms leading to death in the second infection. Overall, 68.8% (11/16) had a similar burden; 18.8% (3/16) had worse symptoms, and 12.5% (2/16) had milder symptoms with the second episode. The conclusion is that, in general, re-infection with different strains is possible and, in some cases, may develop more severe infections with the second episode [44].

Recently, John Ashton, in his review, claimed that the dominance of politics over science as manifested by the rush to abandon all measures of virus control has led to the emergence of a dominant domestic narrative that the pandemic is over. The assumption has been spreading that COVID-19 will only be of nuisance value on a par with the flu in the future. Such a decision has already been taken in Spain. In addition, this theory plays down the associated morbidity and mortality associated with COVID-19 compared to influenza [45].

A similar reaction came from EMA’s vaccine strategy chief Marco Cavaleri almost simultaneously. They also discussed what would be better to do in such a situation. The question is: How should we be thinking about “herd immunity” to COVID-19? According to today’s data for the rate of “herd immunity”, we can not eliminate SARS-CoV-2, and the virus continues to circulate. We do not have enough effective vaccines, and due to the high mutation rate of the virus, the vast majority of the population still can be exposed. That is why we should keep working, predicting and monitoring the situation according to this reality!

Once again, we note that all results in this study were obtained using official data sets. A debatable question is what is the true number of infections with SARS-CoV-2. We cannot answer this question, but we suppose that the most accurate among the reported data are the numbers of COVID-19 deaths. Here we will discuss the results of our study on mortality from COVID-19 and additional mortality during the ongoing pandemic in Bulgaria. Using the available data for deaths in Bulgaria (see [35]) for the last five years 2015, 2016, 2017, 2018 and 2019 before the beginning of the COVID-19 pandemic, we calculate the expected and excess deaths per 100,000 population for the considered period of the pandemic.

An interesting observation is that the times of the excess deaths peaks and the peaks of COVID-19 deaths in Bulgaria coincide (see Figure 14). We will note that the excess deaths during the peak periods are more than twice as big as COVID-19 deaths. This could mean that the COVID-19 is about twice as deadly or that their measures taken to limit its spread cause death as much as the virus itself. In the same time, the WHO reported (see [46]) that the deaths associated with the COVID-19 pandemic between 1 January 2020 and 31 December 2021 were approximately 14.9 million. Hence, globally, the excess deaths from COVID-19 pandemic is about three times more than the reported. Further study is warranted to investigate how the different restriction measures and the vaccination strategy used affect the transmission rate and COVID-19 mortality rate in Bulgaria or in other countries with high mortality.

9. Conclusions

In the present paper, we formulate a new time-dependent deterministic model SIRS-VB with vaccination and vital dynamics. Since available measurements were made at the end of each day of the pandemic, we introduce a semi-implicit finite difference scheme with a step size of 1 day, and then we provide an algorithm for identification of the parameters in the obtained discrete problem. Furthermore, we conduct numerical experiments which show that the calculated parameters values are very close to the values of parameters in the original differential problem. This allows us to calculate realistic values for all compartments in the considered deterministic model by using the suggested identification procedure. The presented study cannot predict the emergence of a new variant or strain of the virus or a new wave caused by it. However, the results in Section 7 show that, when a new wave appears, the model can be efficiently used for making 14-day forecasts for daily numbers of the new infectious cases and the new deaths. The presented prediction methodology is based on the cases/weekdays dependence of the real-time existing data. It requires a very careful analysis of the available data, taking into account the atypical daily values of the parameters during recent weeks and whether such type of values are expected during the forecast period.