A Study on COVID-19 Incidence in Europe through Two SEIR Epidemic Models Which Consider Mixed Contagions from Asymptomatic and Symptomatic Individuals

Abstract

The impact of the SARS-CoV-2 (COVID-19) on the world has been partially controlled through different measures of social isolation and prophylaxis. Two new SEIR (Susceptible-Exposed-Infected-Recovered) models are proposed in order to describe this spread through different countries of Europe. In both models the infectivity of the asymptomatic period during the exposed stage of the disease will be taken into account. The different transmission rates of the SEIR models are calculated by considering the different locations and, more importantly, the lockdown measures implemented in each region. A new classification of these intervention measures will be set and their influence on the values of the transmission rates will be estimated through regression analysis.

1. Introduction

Although prediction of outbreaks such as the one provoked by SARS-CoV-2 (COVID-19), declared a Public Health Emergency of International Concern by the WHO on early 2020, had been made numerous times [1,2], the unexpected spread of this one affected the lives of millions of people as it spread around the world. The rise of diagnosis requiring intensive care and hospitalization had a shocking impact in the healthcare infrastructure in many countries as the numbers went up [3,4,5,6,7]. In March the outbreak, spanning 112 countries and regions, was declared a Pandemic by WHO, who recommended taking action to change the course of the epidemic. The prevention of the propagation of the virus required a method for describing and predicting the apparition of newly infected individuals. Thus, while the main research focused on finding better treatments to the illness and/or boosting the immunity of the individuals through new vaccines, the general mechanism of contagion should also be understood quantitatively, in order to help in the decisions over the taken actions [8,9,10]. Many mathematical models have been proposed in order to understand the dynamic associated to the affected population and their interaction with the virus [11,12,13,14,15,16,17,18,19,20]. There are many different compartmental models that we can apply to the current epidemic [21,22]. The individuals are separated into subpopulations depending on their state regarding the infection, such as the SI (susceptible-infected) model [23], the SIR model (susceptible-infected-recovered) [24], or more complex models when considering more sophisticated reactions to the disease [25]. In this case, the special properties in the propagation of this disease, in which the large incubation (infectious and asymptomatic) time is fundamental to the spread of the disease, can be adjusted to a SEIR model, with two different groups of infected individuals: asymptomatic exposed (E) and the symptomatic infectious (I) apart from the susceptible (S) and recovered (R). There is data being gathered at the time of writing this paper to estimate the number of exposed individuals from the infectious population, as the incubation and recovery time have been properly studied. The parameters associated to the new SEIR model and the the dynamics of social interactions at the affected subpopulations will be then the main objects of study. The transmission mechanism and the role played by asymptomatic individuals in the spread of the disease imply that social-distancing control measures applied in society are effective to fight against the epidemic. In addition, the popularization of diverse non-pharmaceutical interventions (NPI) programs have greatly contributed over time to the awareness of infection-risk behavior which can be possibly averted [26,27]. This subject as well as the different levels of self isolation are studied in this paper. Regrettably, the methods for gathering the data over the different affected countries and regions are not standardized, and the testing capabilities of COVID-19 have not increased at the same rate over different regions and countries. The data describing the impact of the epidemic has shown a great variation depending on the different methods of diagnostics, treatment and surveillance of the affected population [28]. Even the politics may play a part in the transparency of the provided data [29,30,31]. We have chosen to focus our study in regions of the European Union, as we consider that they share similar health infrastructures compared to other parts of the world [32,33], and thus, their vigilance of the epidemic spread will be similar. The set of European countries where the available information is consistent enough to allow us to extract the data required to use in a SEIR model is composed of Spain [34], Italy [35], France [36], the UK [37], Germany [38], Portugal [39], Norway [40] and Sweden [41]. We have created a new classification method regarding the way the NPI are implemented in different regions, so that all countries observe the same standardized levels of isolation and NPI. Thus, all countries can be compared in order to find common patterns on the impact of the control strategies. The network of interactions that define the spread of an infectious disease is usually modeled with differential and discrete equations [42,43]. It involves different types of infected subpopulations, or susceptible to be infected [44,45,46,47], as well as the transitions between them. Although, traditionally, the biggest determinant of the virality of a disease is the Reproduction Number [45,48], which in this case summarizes the impact of the different infectivity rates of the disease, we have chosen to separate these parameters and study them individually. The Reproduction Number does an excellent job at describing the advance of the disease in a population but does not give a direct insight of the causes of contagion (i.e., the contact rate between the different individuals). This will give us the opportunity to focus on the parts that can be directly controlled by the NPI and emphasize the impact of the interaction between individuals on the prevention of the disease. Two different discrete SEIR models, each one describing a different transition from the exposed to the infectious subpopulation, will be introduced. The data of regions of Europe at different organization entity levels (Provinces, Regions or Autonomous Communities and Countries as a whole) will be processed, taking into account the different time intervals depending on the NPI-stages in which they are classified. While the possibility of different strains of the virus spreading at the same time is high, the general assumption will be that most differences regarding the average time of incubation of the virus or the average recovery of an infectious individual will be negligible. Then, the only controllable parameter through NPI control measures and comparable between different regions and stages of social distancing will be the infectivity rate, and it will be calculated given these two models. This paper is structured as follows: in the first two sections, we have performed this brief introduction and present the SEIR models which we will use during our analysis. In Section 3 we will further study theoretical boundaries required to guarantee the model used is adequate to describe correctly the disease spreading through the population. Section 4 will discuss the available data and the extraction of the information necessary to obtain the parameters of the transmission of the disease at the selected intervals and regions. Finally, the results and conclusions derived from them will be presented in Section 5 and Section 6.

2. The $\mu$-SEIR Model

We will explore the dynamics of a compartmental model with four classes to separate the population, corresponding to the subpopulation susceptible to be infected (S) and the infected, which is composed of the subpopulation of individuals which have been exposed (E) with the disease but have not shown the symptoms yet, the symptomatic subpopulation of infected individuals (I) diagnosed with the disease and/or showing symptoms and the recovered subpopulation (R), which presents an immunity which can or cannot be temporary and where individuals are presumed fully immune. The dynamics of the subpopulations can be visualized in the Figure 1:

The model dynamics may be rewritten with the following equations:

$$\begin{array}{c}\hfill S_{i+1}=(1-{\beta }_I{I}_i-{\beta }_E{E}_i)S_i\end{array}$$

(1)

$$\begin{array}{c}\hfill E_{i+1}=({\beta }_I{I}_i+{\beta }_E{E}_i)S_i+(1-\mu )E_i\end{array}$$

(2)

$$\begin{array}{c}\hfill I_{i+1}=\mu E_i+(1-\gamma )I_i\end{array}$$

(3)

$$\begin{array}{c}\hfill R_{i+1}=\gamma I_i+R_i\end{array}$$

(4)

for any integer $i\in {\mathbb{Z}}_{0^{+}}={\mathbb{Z}}_{+}\cup \left\{0\right\}$ and any given initial non-negative conditions $S_0\ge 0$, $E_0\ge 0$, $I_0\ge 0$ and $R_0\ge 0$. The infectivity rates ${\beta }_I$ and ${\beta }_E$ are defined through the force of infection ${\lambda }_i$, which is the rate at which susceptible individuals contract the infection per capita. The value of new infected individuals would be ${\lambda }_i{S}_i$, being $S_i$ the susceptible subpopulation at that moment. The number of infectious individuals will directly influence the force, so the infectivity rates ${\beta }_I$ and ${\beta }_E$ will be introduced in the definition of the force of infection as ${\lambda }_i={\beta }_I{I}_i+{\beta }_E{E}_i$. These rates will depend on multiple factors, such as the average number of contacts any individual of the susceptible subpopulation encounters during certain time with individuals of the exposed and infectious subpopulations, respectively, and the probability of transmission of the disease in a contact between a susceptible individual and an infectious or exposed one [45]. The parameter $\mu$ will represent the average transition rate from the a non-asymptomatic infected subpopulation to the diagnosed/symptomatic one. It is assumed that all exposed individuals recover after experimenting observable symptoms, so the transition to recovered comes uniquely from the infected subpopulation at an average rate $\gamma$. The model will be discrete as the information regarding the spread of the disease is taken and published periodically. Note that (a) the rationale of the sampling period interpretation is that it is unity, typically one day or one week for a correct practical use of the model. (b) the exposed along the incubation period are asymptomatic, including those who then become symptomatic. (c) both the exposed and the infectious produce contagions on the Susceptible which reflects what, in fact, happens with the COVID-19 disease transmission mechanism. The model parameters should be expressed in values of dimensionality being the inverse of the sampling period units. We have chosen the daily period because the available data in the countries we intend to study is released at that rate. It is argued that this model may be of interest for the evaluation of the COVID-19 disease propagation in its blowing-up phase because of the following reasons:

• It is of simple structure and of a discrete nature. Furthermore, it does not need the incorporation of a discretization modulation parameter to guarantee the non-negativity of the sequence solution as other epidemic discrete models usually need. See, for instance, References [24,25] and some references therein.
• It has two potentially distinct coefficient transmission rates (${\beta }_E$ and ${\beta }_I$) which allow to consider infectivity of the susceptible subpopulation from both infective subpopulations: the exposed and infectious subpopulations. This might be potentially advantageous for its use for description of COVID-19 since it is now known that this disease has an infective period at the end of the incubation period and another one along the first days of the symptomatic infectious period.
• It is not very relevant for the current studies related to COVID-19 to evaluate the possible existence of an endemic steady state since the disease is now spreading in very fast phases, to which follows the different social distancing and prophylactic measures taking place in most countries.
• It is later seen from the mathematical study of this model that the susceptible subpopulation is a decreasing sequence and that the exposed and infectious subpopulations increase for consecutive samples under certain conditions of the disease parameters and upper bounds for the susceptible subpopulation. This behavior is also of interest for the use of the model to describe the infection evolution as the disease blows up along a transient period of time.

Boundaries

The subsequent two results are direct from the combination of Equations (1)–(4) and they are respectively concerned with the fact that the total population is constant and the infection spreads, strictly increasing along consecutive samples, under certain constraints of the parameters and the susceptible subpopulation:

Proposition 1.

The total population $N_i=S_i+E_i+I_i+R_i$ is constant for any integer $i\in {\mathbf{Z}}_{{\mathbf{0}}^{+}}$.

Proof. 

When summing up Equations (1)–(4), one gets directly that

$$N_{i+1}=N_i=N_0$$

for any integer $i\in {\mathbf{Z}}_{{\mathbf{0}}^{+}}$. □

The subsequent result gives explicit conditions for the infection to expand at each two consecutive samples provided that the susceptible exceed a certain minimum threshold. This fact intuitively agrees with the known property of the significant verified pandemic spreading in its initial phase when the susceptible levels were very high.

Proposition 2.

Assume that $S_i>max\{\frac{\mu }{{\beta }_E},\frac{\gamma (1-\gamma )}{\mu {\beta }_I}\}$. Then $E_{i+1}>E_i$ and $I_{i+1}>I_i$$\forall i>1$.

Proof. 

Define $\Delta E_i=E_{i+1}-E_i$ and $\Delta I_i=I_{i+1}-I_i$. From Equation (2) we obtain that $\Delta E_i={\beta }_I{I}_i{S}_i+({\beta }_E{S}_i-\mu )E_i\ge {\beta }_I{I}_i{S}_i+({\beta }_E{S}_i-\mu )E_i$ which, being $S_i>\frac{\mu }{{\beta }_E}$ implies that $\Delta E_i>0$ and $E_{i+1}>E_i$. Now, from (3) we obtain that $\Delta I_i=\mu E_i-\gamma I_i$. Then $\Delta I_{i+1}=\mu E_{i+1}-\gamma I_{i+1}=\mu (E_i+\Delta E_i)-\gamma (I_i+\Delta I_i)=(1-\gamma )\Delta I_i+\mu \Delta E_i$ Which would be positive for all i+1 if $\Delta E_i>0$, which it has already be proven, and $\Delta I_i\ge 0$. If $\Delta I_i<0$, $\Delta I_{i+1}>0$ would be true only if $\Delta E_i>\frac{\gamma -1}{\mu }\Delta I_i$. From Equations (2) and (3) we get that

$$\begin{array}{c}\hfill {\beta }_I{S}_i{I}_i+E_i({\beta }_E{S}_i-\mu )>\frac{\gamma -1}{\mu }(\mu E_i-\gamma I_i)\end{array}$$

(5)

$$\begin{array}{c}\hfill E_i({\beta }_E{S}_i+1-\gamma -\mu )>I_i\left(\frac{\gamma (1-\gamma )}{\mu }-{\beta }_I{S}_i\right)\end{array}$$

(6)

Given that $E_i>0,I_i>0$ $\forall i\in {\mathbf{Z}}_{\mathbf{0}+}$, these inequalities will be always true if

$$\left\{\begin{array}{ccc}{S}_i\hfill & >\hfill & \frac{\gamma +\mu -1}{{\beta }_E}\hfill \\ S_i\hfill & >\hfill & \frac{\gamma (1-\gamma )}{\mu {\beta }_I}\hfill \end{array}\right.$$

(7)

Then, if $S_i>max\{\frac{\gamma (1-\gamma )}{\mu {\beta }_I},\frac{\mu +\gamma -1}{{\beta }_E},\frac{\mu }{{\beta }_E}\}=max\{\frac{\gamma (1-\gamma )}{\mu {\beta }_I},\frac{\mu }{{\beta }_E}\}$ the proposition holds. □

The vector sequence solution of dimension four is non-negative for all samples for any given non-negative initial conditions of all the subpopulations under some reasonable extra constraints on the initial conditions and parameters (so as to keep jointly both stability and non-negativity of the solution) as proved in the subsequent result:

Theorem 1.

Assume that $\mu ,\gamma \in (0,1]$, and ${\beta }_I/{\beta }_E\le \gamma /\mu$ with ${\beta }_I>0$ and ${\beta }_E>0$. Assume also that the initial conditions $(S_0,E_0,I_0,R_0)$ are subject to

$$\begin{array}{c}\hfill \left\{\begin{array}{ccc}0\hfill & \le \hfill & S_0\le \frac{\mu }{2{\beta }_E}\hfill \\ 0\hfill & \le \hfill & E_0\le \frac{1}{2{\beta }_E}\hfill \\ 0\hfill & \le \hfill & I_0\le \frac{1}{2{\beta }_I}\hfill \\ 0\hfill & \le \hfill & R_0\hfill \end{array}\right.\end{array}$$

(8)

Then, the vector sequence solution is non-negative for all samples. In addition, the sequences of all the subpopulations and that of the total population are bounded.

Proof. 

The proof follows by complete induction. Assume that there exists an integer $k\ge 0$ such that $min(S_i,E_i,I_i,R_i)\ge 0$, $0\le S_i\le \frac{\mu }{2{\beta }_E}$, $0\le E_i\le \frac{1}{2{\beta }_E}$ and $0\le I_i\le \frac{1}{2{\beta }_I}$ for any $i=0,1,2,\dots ,k$. Such an integer $k\ge 0$ always exists since above constraints hold from the given hypotheses on the initial conditions (at least), i.e., for $k=0$. It has to be proved that the constraints still hold for $k+1$. First note that, since $E_k\le \frac{1}{2{\beta }_E}$, $I_k\le \frac{1}{2{\beta }_I}$ and $0\le S_k$, one has:

$$S_{k+1}=(1-{\beta }_I{I}_k-{\beta }_E{E}_k)S_k\ge 0$$

(9)

and, since $E_k\ge 0$ and $I_k\ge 0$, one has that

$$I_{k+1}=(1-\gamma )I_k+\mu E_k\ge 0$$

(10)

On the other hand, since $E_k\le \frac{1}{2{\beta }_E}$, $I_k\le \frac{1}{2{\beta }_I}$, $0\le S_k\le \frac{\mu }{2{\beta }_E}$, $\gamma \in (0,1]$ and ${\beta }_I/{\beta }_E\le \gamma /\mu$, one has that:

$$\begin{array}{cc}\hfill 0\le E_{k+1}& =(1-\mu )E_k+({\beta }_I{I}_k+{\beta }_E{E}_k)S_k\hfill \\ \hfill & \le \frac{1-\mu }{2{\beta }_E}+S_k\le \frac{1-\mu }{2{\beta }_E}+\frac{\mu }{2{\beta }_E}\hfill & \\ \hfill & =\frac{1}{2{\beta }_E}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill 0\le I_{k+1}& =(1-\gamma )I_k+\mu E_k\hfill \\ \hfill & \le \frac{1-\gamma }{2{\beta }_I}+\frac{\mu }{2{\beta }_E}\le \frac{1-\gamma }{2{\beta }_I}+\frac{\gamma }{2{\beta }_I}\hfill \\ \hfill & =\frac{1}{2{\beta }_I}\hfill \end{array}$$

In addition, since $S_k\le \frac{\mu }{2{\beta }_E}$ and $0\le (1-{\beta }_I{I}_k-{\beta }_E{E}_k)$ since $0\le E_k\le \frac{1}{2{\beta }_E}$; $0\le I_k\le \frac{1}{2{\beta }_I}$, it follows that

$$0\le S_{k+1}=(1-{\beta }_I{I}_k-{\beta }_E{E}_k)S_k\le S_k\le \frac{\mu }{2{\beta }_E}$$

Finally, since $R_k\ge 0,I_k\ge 0$ and $\gamma \ge 0$ then $R_{k+1}\ge 0$. It has been proved that $(min(S_i,E_i,I_i,R_i)\ge 0;i=0,1,2,\dots ,k)$ implies $(min(S_i,E_i,I_i,R_i)\ge 0;i=0,1,2,\dots ,k+1)$, and $\left(S_i\le \frac{\mu }{2{\beta }_E},E_i\le \frac{1}{2{\beta }_E},I_i\le \frac{1}{2{\beta }_I};i=0,1,\dots ,k\right)$ implies ($S_i\le \frac{\mu }{2{\beta }_E},E_i\le \frac{1}{2{\beta }_E},I_i\le \frac{1}{2{\beta }_I}$; $i=0,1,\dots ,k+1)$. As a result the proof is complete since $min(S_0,I_0,E_0,R_0)\ge 0$ implies $min(S_i,E_i,I_i,R_i)\ge 0;\forall i\in {\mathbf{Z}}_{+}$ and $\left(S_0\le \frac{\mu }{2{\beta }_E},E_0\le \frac{1}{2{\beta }_E},I_0\le \frac{1}{2{\beta }_I}\right)\to S_i\le \frac{\mu }{2{\beta }_E}$, $E_i\le \frac{1}{2{\beta }_E},I_i\le \frac{1}{2{\beta }_I};\forall i\in {\mathbf{Z}}_{+}$ which implies also that the sequences of all the subpopulations and that of the total population are bounded. □

Remark 1.

The given constraints on the upper bounds of the subpopulations which guarantee that Theorem 1 holds are reasonable since the disease transmission rates are typically small in normalized models and much smaller in un-normalized ones. Typically, the normalized disease transmission rate is normalized by the total population in true-mass action epidemic models. On the other hand, the constraint ${\beta }_I/{\beta }_E\le \gamma /\mu$ is reasonable, since typically $\gamma <\mu$ and the transmission rate of the infectious individuals are on average less than that of the exposed ones since the strongest infective period happens at the end period of the disease incubation. In addition, note that if ${\beta }_I={\beta }_E=\beta$ then the constraints $E_i\le \frac{1}{2{\beta }_E}$, $I_i\le \frac{1}{2{\beta }_I}$ of Theorem 1 become a single constraint $E_i+I_i\le \frac{1}{\beta };\forall i\in {\mathbf{Z}}_{\mathbf{0}+}$ and Theorem 1 remains valid. It also remains valid if $max(E_i,I_i)<\frac{1}{{\beta }_I+{\beta }_E};\forall i\in {\mathbf{Z}}_{\mathbf{0}+}$.

The subsequent technical result is useful later on:

Proposition 3.

Assume that the hypotheses of Theorem 1 hold. Thus, if $\underset{i\to \infty }{lim}{E}_i=0$ then ${\displaystyle \underset{i\to \infty }{lim}{I}_i=0}$ and, conversely, if ${\displaystyle \underset{i\to \infty }{lim}{I}_i=0}$ then $\underset{i\to \infty }{lim}{E}_i=0$.

Proof. 

One gets from Equation (3) via recursive calculations that

$$I_{i+1}={(1-\gamma )}^{i+1}{I}_0+\mu \sum _{j=0}^i{(1-\gamma )}^{i-j}{E}_j$$

so that $\underset{i\to \infty }{lim}{E}_i=0$ implies that for any given arbitrary real constant $ϵ>0$, there exists a non-negative integer $m_0=m\left(ϵ\right)$ such that $\underset{\infty >j\ge m_0}{sup}{E}_j\le ϵ$ and $lim\underset{i\to \infty }{sup}\left(\underset{i-m_0\le j\le i}{sup}{E}_j\right)\le ϵ$ then one has for any integer $m\ge m_0$ that

$$\begin{array}{ccc}\hfill I_{i+1}& =& {(1-\gamma )}^{m+1}{I}_{i-m}+\mu \sum _{j=i-m}^i{(1-\gamma )}^{i-j}{E}_j\hfill \\ & \le & {(1-\gamma )}^{m+1}{I}_{i-m}+\mu ϵ\sum _{j=i-m}^i{(1-\gamma )}^{i-j}\hfill \\ & \le & {(1-\gamma )}^{m+1}{I}_{i-m}+\frac{\mu ϵ}{\gamma }\hfill \end{array}$$

so that $\underset{(i-m)\to \infty }{\underset{\underset{i\to \infty }{m\to \infty }}{lim}}{I}_{i+1}=\mu ϵ/\gamma$ for any integer $m\ge m_0$. Since $ϵ>0$ may be chosen arbitrarily small then $m_0=\underset{ϵ\to 0}{lim}{m}_0\left(ϵ\right)=\infty$ and so $m(\ge m_0)\to \infty$, then ${\displaystyle \underset{i\to \infty }{lim}{I}_i=0}$. The first implication has been proved. It remains to be proved that, if $\underset{i\to \infty }{lim}{I}_i=0$, then $\underset{i\to \infty }{lim}{E}_i=0$. One has from Equation (2) that:

$$E_{i+1}=(1-\mu +{\beta }_E{S}_i)E_i+{\beta }_I{S}_i{I}_i$$

Since $\underset{i\to \infty }{lim}{I}_i=0$ and ${\left\{S_i\right\}}_{i=0}^{\infty }$ is a bounded sequence then $\underset{i\to \infty }{lim}{\beta }_I{I}_i{S}_i=0$, and then $\underset{i\to \infty }{lim}(E_{i+1}-(1-\mu +{\beta }_E{S}_i)E_i)=0$. Two cases can arise, namely:

1.

$\underset{i\to \infty }{lim}{E}_i=0$ and the proof finishes,

2.

$\underset{i\to \infty }{lim}{E}_i\ne 0$. Then $lim\underset{i\to \infty }{inf}{E}_i>0$ so that, from Equation (1), and since ${\displaystyle \underset{i\to \infty }{lim}{I}_i=0}$ and $\frac{1}{2{\beta }_E}\ge E_i\ge 0\forall i\ge 0$ from Theorem 1, one has that

$$lim\underset{i\to \infty }{sup}\left(1-{\beta }_I{I}_i-{\beta }_E{E}_i\right)=lim\underset{i\to \infty }{sup}(1-{\beta }_E{E}_i)=\left(1-{\beta }_{E}lim\underset{i\to \infty }{inf}{E}_i\right)=a$$

with $a\in [\frac{1}{2},1)$, which implies from Equation (1) and Theorem 1 that $0=lim\underset{i\to \infty }{sup}\left(S_{i+1}-a^{i+1}{S}_0\right)=\underset{i\to \infty }{lim\; sup}{S}_{i+1}$ so that there exist

$$\underset{i\to \infty }{lim}{S}_{i+1}=lim\underset{i\to \infty }{sup}{S}_{i+1}=0$$

Then, using again $E_{i+1}=\left(1-\mu +{\beta }_E{S}_i\right)E_i+{\beta }_I{S}_i{I}_i$, with $\underset{i\to \infty }{lim}{I}_i=0$ and $\underset{i\to \infty }{lim}{S}_i=0$ yields $\underset{i\to \infty }{lim}\left(E_{i+1}-(1-\mu )E_i\right)=0$ which, in turn, implies that:

(a)

Either $\underset{i\to \infty }{lim}{E}_i=0$ and the result is proved,

(b)

or $\underset{i\to \infty }{lim}{E}_i=E^{*}\ne 0$ which implies that $1=1-\mu$, so $\mu =0$, which contradicts $\mu \in (0,1]$ and thus, this case cannot occur

(c)

or $\frac{E_{i+1}}{E_i}\to 1-\mu \in [0,1)$. Then, there is a subsequence ${\left\{E_i\right\}}_{i=N}^{\infty }$ of ${\left\{E_i\right\}}_{i=0}^{\infty }$, for some finite integer N, which is non-negative and strictly decreasing so that $\underset{i\to \infty }{lim}{E}_i=0$ and the result is proved.

The proposition has been fully proved. □

By inspection of Equations (1)–(4), it follows the existence of a, in general, non-unique disease-free equilibrium point. Thus,

Proposition 4.

There is a disease-free equilibrium point $x_{dfe}={(S_{dfe},0,0,R_{dfe}=N_0-S_{dfe})}^T$ with $S_{dfe}$, and then $N_{dfe}$, being dependent on the initial conditions.

The demonstration is trivial from Equations (1)–(4) by writing $E_{i+1}=E_i=I_{i+1}=I_i=0=E_{dfe}=Idfe$ and $N_{dfe}=S_{dfe}+R_{dfe}$ The following result is concerned with the convergence of the vector sequence solution to the disease-free equilibrium point if the transmission rate of the exposed (E) exceeds that of the infected (I) under a certain minimum lower-bound constraint. This is a reasonable consideration since the exposed subpopulation in the model is also infective and the infectivity along the incubation period is normally of the same magnitude of the infectious subpopulation [49,50]. Basically, if the constraint ${\beta }_E/{\beta }_I\ge \mu /\gamma$ of Theorem 1 for $\gamma \in (0,1)$, then the properties of non-negativity of the solution and its asymptotic convergence to the disease-free equilibrium point are jointly guaranteed under the constraints for the initial conditions of Theorem 1.

Theorem 2.

Assume that $\mu \in (0,1],\gamma \in (0,1]$, ${\beta }_I>0$ and ${\beta }_E>0$ and that conditions of Theorem 1 hold. Assume also that the initial conditions are subject to $0\le S_0\le \frac{\mu }{2{\beta }_E}\le \frac{\gamma \mu }{\mu {\beta }_I+\gamma {\beta }_E},0\le E_0\le \frac{1}{2{\beta }_E},0\le I_0\le \frac{1}{2{\beta }_I},\mathrm{and}0\le R_0$. Then, the vector sequence solution is non-negative and bounded for all samples and converges asymptotically to the disease-free equilibrium at exponential rate.

Proof. 

Equations (2) and (3) can jointly be expressed as follows:

$$\left[\begin{array}{c}{E}_{i+1}\\ I_{i+1}\end{array}\right]=\left[\begin{array}{cc}1-{\beta }_E{S}_i-\mu & {\beta }_I{S}_i\\ \mu & 1-\gamma \end{array}\right]\left[\begin{array}{c}{E}_i\\ I_i\end{array}\right]$$

(11)

The characteristic equation of the discrete system from Equation (11) for the case when $S_i$ is constant, that is, assuming that the system (11) is time-invariant is

$$\begin{array}{cc}\hfill q_i\left(z\right)=& q_i(z,S_i)=(z-1+\mu +{\beta }_E{S}_i)(z+\gamma -1)-\mu {\beta }_I{S}_i\\ \hfill =& z^2+(\mu +\gamma +{\beta }_E{S}_i-2)z+(\gamma -1)(\mu -1+{\beta }_E{S}_i)-\mu {\beta }_I{S}_i\end{array}$$

(12)

The system Equation (11), provided invariant, is stable if its modes are in $\left|z\right|<1$. The conditions can be discussed by using the bilinear transformation $z=\frac{1+s}{1-s}$ which is a comfortable transformation (except for $s=1$, $z=+\infty$) which transforms the complex region $\left|z\right|\le 1$ into the continuous complex region $Re\left(s\right)\le 0$. This transformation allows one to apply comformable the Routh–Hurwitz stability criterion to the transformed characteristic equation which is simpler to apply, especially for this simple second-order system, than discrete Jury’s stability criterion. Thus, one defines ${\widehat{q}}_{si}\left(s\right)=q_i\left(z=\frac{1+s}{1-s}\right)$ so that its zeros being the solution of ${\widehat{q}}_{si}\left(s\right)=0$ are equivalent, after reducing to the common denominator, to the zeros of

$$\begin{array}{ccc}\hfill q_{si}\left(s\right)& =& {(1+s)}^2+(\mu +\gamma +{\beta }_E{S}_i-2)(1-s^2)\hfill \\ & & +\left[(\gamma -1)(\mu -1)-S_i\left(\mu {\beta }_I-{\beta }_E(\gamma -1)\right)\right]{(1-s)}^2\hfill \\ & =& a\left(S_i\right)s^2+b\left(S_i\right)s+c\left(S_i\right)\hfill \end{array}$$

(13)

where

$$\begin{array}{c}\hfill \begin{array}{ccc}{a}_i=a\left(S_i\right)\hfill & =\hfill & (2-\gamma )(2-\mu )-S_i\left((2-\gamma ){\beta }_E+\mu {\beta }_I\right)\hfill \\ b_i=b\left(S_i\right)\hfill & =\hfill & 2\left[\gamma +\mu -\gamma \mu +S_i(\mu {\beta }_I+(1-\gamma ){\beta }_E\right)]\hfill \\ c_i=c\left(S_i\right)\hfill & =\hfill & \gamma \mu +S_i(\gamma {\beta }_E-\mu {\beta }_I)\hfill \end{array}\end{array}$$

(14)

It is now seen that the above three coefficients are positive for any $i\in {\mathbf{Z}}_{\mathbf{0}+}$ provided that the constraints of Theorem 1 hold, that is, (8), ${\beta }_I/{\beta }_E\le \gamma /\mu$, with $\mu \in (0,1]$, and $\gamma \in (0,1]$. Note that it is guaranteed from Equations (14) that:

1.

$c_i\ge 0;\forall i\in {\mathbf{Z}}_{\mathbf{0}+}$ if $\gamma {\beta }_E-\mu {\beta }_I\ge 0$. Since the non-negativity conditions from Theorem 1 requires that ${\beta }_I/{\beta }_E\le \gamma /\mu$ it follows that $c_i>0$.

2.

$b_i>0$ always, as the conditions $\gamma \in (0,1)$ and $\mu \in (0,1)$ implies that $b_i=2[\gamma +\mu (1-\gamma )+S_i\mu {\beta }_I+S_i(1-\gamma ){\beta }_E]\ge 0$, as it is a sum of positive terms.

3.

$a_i>0;\forall i\in {\mathbf{Z}}_{\mathbf{0}+}$ iff $\frac{(2-\gamma )(2-\mu )}{(2-\gamma ){\beta }_E+\mu {\beta }_I}>S_i$. Since the condition of non-negativity from Theorem 1$\frac{{\beta }_I}{{\beta }_E}\le \frac{\gamma }{\mu }$, or equivalently, ${\beta }_E/{\beta }_I\ge \mu /\gamma$, has to hold then

$$\frac{(2-\gamma )(2-\mu )}{(2-\gamma ){\beta }_E+\mu {\beta }_I}\ge \frac{(2-\gamma )(2-\mu )}{(2-\gamma ){\beta }_E+\gamma {\beta }_E}=\frac{(2-\gamma )(2-\mu )}{2{\beta }_E}>\frac{\mu }{2{\beta }_E}\ge S_i$$

so that $a_i>0$.

It has been proved that the three coefficients of equations in (14) are positive for any $i\in {\mathbb{Z}}_{0+}$. This implies that all the polynomials of $q_{si}\left(s\right)$ Equation (13) satisfy the Hurwitz’s test for any $i\in {\mathbf{Z}}_{\mathbf{0}+}$. Since all those polynomials are of second order then all of them are from Routh–Hurwitz criterion (i.e., with both zeros in the open left-hand-side complex plane $Res<0$). Then, all the discrete polynomials $q_i\left(z\right)$ of Equation (12) for any $i\in {\mathbf{Z}}_{\mathbf{0}+}$ are stable (i.e., with both zeros in the open unit circle centered at zero, that is $\left|z\right|<1$). As a result, all the time-invariant second-order discrete systems of Equation (11) (each for an assumed constant $S_i$), for any $i\in {\mathbf{Z}}_{\mathbf{0}+}$ are exponentially stable. Thus, each norm of matrix of dynamics:

$${\mathbf{A}}_i=\mathbf{A}\left(S_i\right)=\left[\begin{array}{cc}1+{\beta }_E{S}_i-\mu & {\beta }_I{S}_i\\ \mu & 1-\gamma \end{array}\right]$$

(15)

whose eigenvalues are the zeros of their corresponding polynomial $q_i\left(z\right)$ for each $i\in {\mathbb{Z}}_{0+}$, satisfies $\left|\right|{\mathbf{A}}_i{\left|\right|}^k\le K_i{\rho }_i^k$ for some real constant $\infty >K_i\ge 1$, which is norm dependent, and ${\rho }_i\in (0,1)$ which is the spectral radius of ${\mathbf{A}}_i$: That is the maximum absolute value of its two eigenvalues. Now define the real numbers:

$$K=lim\underset{n\to \infty }{sup}\left(\underset{0\le i\le n}{sup}{K}_i\right)\in [1,\infty );\rho =lim\underset{n\to \infty }{sup}\left(\underset{0\le i\le n}{sup}{\rho }_i\right)\in (0,1)$$

(16)

As a result, one gets from Equation (11) that:

$$\begin{array}{cc}\hfill ∥\begin{array}{c}{E}_i\\ I_i\end{array}∥\le & \parallel A_{i-1}\parallel ∥\begin{array}{c}{E}_{i-1}\\ I_{i-1}\end{array}∥\\ \hfill \le & \parallel A_{i-1}\parallel \parallel A_{i-2}\parallel ∥\begin{array}{c}{E}_{i-2}\\ I_{i-2}\end{array}∥\\ \hfill \le & max\left({\parallel A_{i-1}\parallel }^2,{\parallel A_{i-2}\parallel }^2\right)∥\begin{array}{c}{E}_{i-2}\\ I_{i-2}\end{array}∥\\ \hfill \le & {\displaystyle \underset{0\le j\le i-1}{sup}\left({\parallel A_j\parallel }^i\right)∥\left(\begin{array}{c}{E}_0\\ I_0\end{array}\right)∥}\\ \hfill \le & K{\rho }^i∥\left(\begin{array}{c}{E}_0\\ I_0\end{array}\right)∥=0\end{array}$$

(17)

Then, $\underset{i\to \infty }{lim}{E}_i=\underset{i\to \infty }{lim}{I}_i=0$ asymptomatically at exponential rate and the solution sequences ${\left\{E_i\right\}}_{i=0}^{\infty }$, ${\left\{I_i\right\}}_{i=0}^{\infty }$, ${\left\{S_i\right\}}_{i=0}^{\infty }$, ${\left\{R_i\right\}}_{i=0}^{\infty }$ of Equation (1)–(4) are non-negative and bounded from Theorem 1 for any given finite initial conditions which satisfy Equation (8). In addition, the map of $\left[0,\frac{\mu }{2{\beta }_E}\right]$ to itself defined by Equation (1) is non-expansive so that it has a fixed point which is the first component $S_{dfe}$ of the disease-free equilibrium point and $\underset{i\to \infty }{lim}{S}_i=S_{dfe}$. Finally, $\underset{i\to \infty }{lim}{R}_i=R_{dfe}=N_0-S_{dfe}$ since $\underset{i\to \infty }{lim}{E}_i=\underset{i\to \infty }{lim}{I}_i=0$ and the proof is complete. □

Note that the given model is claimed for its usefulness for short-term predictions in the evolution phase when the disease is blowing up. It is neither considered that vaccination is available for use or that there is immunity loss allowing for the increase again of the susceptible numbers after a certain delay. Therefore, the susceptible subpopulation is given by a decreasing sequence. It is now proved that the proposed model does not have an endemic equilibrium point.

Theorem 3.

Assume that $\mu ,\gamma \in (0,1]$. Then the delay-free model Equations (1)–(4) has no endemic equilibrium point $x_{end}={(S_{end},E_{end},I_{end},R_{end}=N_0-(S_{end}+E_{end}+I_{end}))}^T$.

Proof. 

Assume that there is an endemic equilibrium point with $S_{end}>0$. Then, from Equation (1), one gets ${\beta }_I{I}_{end}+{\beta }_E{E}_{end}=0$ so that $E_{end}=I_{end}=0$ which concludes that the eventual equilibrium point is the disease-free one. Thus, in order to the endemic equilibrium to exist, $S_{end}=0$ and ${\beta }_I{I}_{end}+{\beta }_E{E}_{end}=c>0$. From $S_{end}=0$ and Equation (2), one has that $0=S_{end}c-\mu E_{end}=-\mu E_{end}$. Then, from Equation (3), one gets $I_{end}=(1-\gamma )I_{end}+\mu E_{end}=(1-\gamma )I_{end}$, which implies that $I_{end}=0$ if $\gamma \in (0,1)$ and $I_{end}>0$ is arbitrary if $\gamma =1$. As a result, $E_{end}=I_{end}=0$ if $\gamma \in (0,1)$ and the equilibrium point is the disease-free one. Assume then that $\gamma =1$, $E_{end}=0$ and $I_{end}>0$ is arbitrary. Thus, one gets from Equation (4) that $R_{end}=R_{end}+I_{end}$, which implies that $I_{end}=0$, a direct contradiction to $I_{end}>0$.

As a result, it is concluded that the endemic equilibrium point does not exist. □

The subsequent remark discusses the fact that under proportional linear feedback vaccination with a rather weaker condition on the control gain, the disease-free susceptible level can be zeroed so that all the population results asymptotically recovered in the disease-free equilibrium.

Remark 2.

Note that the given model does not guarantee the convergence to zero of the susceptible subpopulation, that is, it does not guarantee that $S_{dfe}=0$ which would imply as a result that $R_{dfe}=N_0$; namely, it is not guaranteed that all the population would become asymptotically recovered. The reason is that the coefficient sequence in parenthesis in Equation (1) converges asymptotically to one when $i\to \infty$ as the infection vanishes so that the solution sequence converges asymptotically to the disease-free equilibrium point. This is not very surprising since former studies indicate that a small residual permanent nonzero susceptible population can exist even if the infection asymptotically extinguishes, implying that not all the population becomes asymptotically recovered. See, for instance, [51,52]. The reason is that the coefficient sequence of Equation (1) of general term $a_i=1-({\beta }_I{I}_i+{\beta }_E{E}_i);i\in {\mathbb{Z}}_{0+}$, leading to $S_{i+1}=a_i{S}_i$ verifies $a_i\in (0,1)$ for $i\in {\mathbb{Z}}_{0+}$ provided that $min(E_0,I_0)>0$ with ${\left\{a_i\right\}}_0^{\infty }\to 1$ since ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$ and ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$. However, if a feedback vaccination law proportional to the susceptible of the form $V_i=k_i{S}_i$ is incorporated to the dynamic Equations (1)–(4), so that $k_i{S}_i$ are removed from the right-hand side of Equation (1) and added to the right-hand side of Equation (4), then the model becomes modified as follows:

$$\begin{array}{ccc}\hfill S_{i+1}& =& (1-(k_i+{\beta }_I{I}_i+{\beta }_E{E}_i))S_i\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill E_{i+1}& =& (1-\mu )E_i+({\beta }_I{I}_i+{\beta }_E{E}_I)S_i\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill I_{i+1}& =& (1-\gamma )I_i+\mu E_i\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill R_{i+1}& =& R_i+\gamma I_i+k_i{S}_i\hfill \end{array}$$

(21)

Now, if the sequence of control gains ${\left\{k_i\right\}}_{i=0}^{\infty }$ is positive and it satisfies furthermore the rather weak constraint ${\displaystyle lim\underset{i\to \infty }{sup}(k_i+{\beta }_I{I}_i+{\beta }_E{E}_i)<1}$ then $S_{dfe}=0$ and $R_{dfe}=N_0$. Note that in order that ${\displaystyle lim\underset{i\to \infty }{sup}(k_i+{\beta }_I{I}_i+{\beta }_E{E}_i)<1}$, it suffices that some prefixed finite

$$N\in {\mathbb{Z}}_{0+},k_i<1-({\beta }_I{I}_i+{\beta }_E{E}_i);\forall i(\ge N)\in {\mathbb{Z}}_{0+}$$

Note also that the last sufficiency-type condition does not require that the control gain converges to a limit but, if it is the case, such limit has to be less than one.

The vaccination law discussed in the above Remark 2 can be even weakened by certain vaccination laws whose gains asymptotically vanish at smaller rates than exponential ones. The subsequent result discusses a more general vaccination law which includes that of Remark 2 and allows the equilibrium susceptible subpopulation to reach a zero value under an asymptotically vanishing vaccination gain at lower rate than exponential:

Theorem 4.

Assume that the following hypotheses hold:

1.

The sequence of control gains ${\left\{k_i\right\}}_{i=0}^{\infty }$ of the linear feedback vaccination law $V_i=k_i{S}_i;i\in {\mathbb{Z}}_{0+}$ is generated with $k_0<1-({\beta }_I{I}_0+{\beta }_E{E}_0)$ and a general term $k_i=\frac{{ϵ}_i}{i^{{\theta }_i}};\forall i\in {\mathbb{Z}}_{+}$ being subject to the following constraints:

• ${\left\{{\theta }_i\right\}}_{i=0}^{\infty }\subset [0,1)$, with ${\theta }_i\in \left(max\left(0,\frac{ln\left(\frac{{ϵ}_i}{1-({\beta }_I{I}_i+{\beta }_E{E}_i)}\right)}{ln\left(i\right)}\right),1\right);\forall i\in {\mathbb{Z}}_{+}$ and
• ${\left\{{ϵ}_i\right\}}_{i=0}^{\infty }\subset (0,1)$, with ${ϵ}_i\in (0,1-({\beta }_I{I}_0+{\beta }_E{E}_0));\forall i\in {\mathbb{Z}}_{0+}$, having a sub-sequence ${\left\{{ϵ}_i\right\}}_{i=N}^{\infty }\subset \left(i^{{\theta }_i-1},1-({\beta }_I{I}_0+{\beta }_E{E}_0)\right)$ for some finite $N\in {\mathbb{Z}}_{0+}$.

2.

all the hypotheses of Theorem 1 hold with the further restrictions for $E_0$ and $I_0$ that $0\le E_0\le (1-\overline{k})/2{\beta }_E$ and $0\le I_0\le (1-\overline{k})/2{\beta }_I$, where $\overline{k}={sup}_{i\in {\mathbb{Z}}_{0+}}{k}_i$.

Then, ${\sum }_{i=0}^{\infty }{S}_i<+\infty$ and $\underset{i\to \infty }{lim}{S}_i=S_{dfe}=0$.

Proof. 

First, note that $k_i=\frac{{ϵ}_i}{i^{{\theta }_i}}<1-({\beta }_I{I}_i+{\beta }_E{E}_i)$, since ${ϵ}_i<i^{{\theta }_i}(1-({\beta }_I{I}_i+{\beta }_E{E}_i))$;$\forall i\in {\mathbb{Z}}_{+}$ (as ${\theta }_i>\frac{ln\left(\frac{{ϵ}_i}{1-({\beta }_I{I}_i+{\beta }_E{E}_i)}\right)}{ln\left(i\right)}$) and, together with $k_0<1-({\beta }_I{I}_0+{\beta }_E{E}_0)$ yields $k_i<(1-({\beta }_I{I}_i+{\beta }_E{E}_i));\forall i\in {\mathbb{Z}}_{0+}$ (see Remark 2) so that one concludes that ${\left\{S_i\right\}}_{i=0}^{\infty }\subset {\mathbb{R}}_{0+}$. Note also that the stronger constraints $0\le E_0\le \frac{1-\overline{k}}{2{\beta }_E}$ and $0\le I_0\le \frac{1-\overline{k}}{2{\beta }_I}$ over the parallel ones of Theorem 1 imply the non-negativity of all the subpopulation sequences by a direct extension of the proof of Theorem 1 with those further constraints. In addition, note that since ${ϵ}_i\in (0,1-({\beta }_I{I}_i+{\beta }_E{E}_i))$ then the subsequence ${\left\{{ϵ}_i\right\}}_{i=N}^{\infty }\subset (i^{{\theta }_i-1},1-({\beta }_I{I}_i+{\beta }_E{E}_i))$ of ${\left\{{ϵ}_i\right\}}_{i=0}^{\infty }$ exists since ${\theta }_i<1$ implies that the sequence ${\left\{i^{{\theta }_i-1}\right\}}_{i=0}^{\infty }$ is strictly decreasing at a smaller convergence rate than exponential leading to $\underset{i\to \infty }{lim}{i}^{{\theta }_i-1}=0$. Furthermore, since $\underset{i\to \infty }{lim}(E_i,I_i)=0$ at exponential rate, then $\underset{i\to \infty }{lim}(1-({\beta }_I{I}_i+{\beta }_E{E}_i))=1$ and then there exists some finite $N\in {\mathbb{Z}}_{0+}$ such that $i^{{\theta }_i-1}<1-({\beta }_I{I}_i+{\beta }_E{E}_i);\forall i(\ge N)\in {\mathbb{Z}}_{0+}$. Thus, the sequence of nonempty real intervals ${\left\{\left(i^{{\theta }_i-1},1-({\beta }_I{I}_i+{\beta }_E{E}_i)\right)\right\}}_{i=N}^{\infty }$ exists, and the claimed subsequence ${\left\{ϵ\right\}}_{i=N}^{\infty }\subset (i^{{\theta }_i-1},1-({\beta }_I{I}_i+{\beta }_E{E}_i))$ exists as well. Now, note that

$$\left(1-\frac{S_{i+1}}{S_i}\right)=1-1+k_i+{\beta }_I{I}_i+{\beta }_E{E}_i=\frac{{ϵ}_i}{i^{{\theta }_i}}+{\beta }_I{I}_i+{\beta }_E{E}_i$$

so that, since ${\left\{{ϵ}_i\right\}}_{i=N}^{\infty }\subset \left(i^{{\theta }_i-1},1\right)$ one has for any $i(\ge N)\in {\mathbb{Z}}_{0+}$ that

$$1-\frac{S_{i+1}}{S_i}-{\beta }_I{I}_i-{\beta }_E{E}_i=\frac{{ϵ}_i}{i^{{\theta }_i}}>\frac{i^{{\theta }_i-1}}{i^{{\theta }_i}}=\frac{1}{i};\forall i(\ge N)\in {\mathbb{Z}}_{0+}$$

which implies that

$$1-\frac{S_{i+1}}{S_i}-{\beta }_I{I}_i-{\beta }_E{E}_i-\frac{1}{i}>0;\forall i(\ge N)\in {\mathbb{Z}}_{0+}$$

which implies, in turn, that

$$i\left(1-\frac{S_{i+1}}{S_i}-{\beta }_I{I}_i-{\beta }_E{E}_i\right)-1>0$$

and, since ${\displaystyle \underset{i\to \infty }{lim}}({\beta }_I{I}_i+{\beta }_E{E}_i)=0$, one gets:

$$lim\underset{i\to \infty }{inf}\left[i\left(1-{\beta }_I{I}_i-{\beta }_E{E}_i-\frac{S_{i+1}}{S_i}\right)\right]-1=lim\underset{i\to \infty }{inf}\left[i\left(1-\frac{S_{i+1}}{S_i}\right)\right]-1>0$$

which is equivalent to ${\displaystyle \underset{i\to \infty }{lim}}inf\left[i\left(1-\frac{S_{i+1}}{S_i}\right)\right]-1>0$ and, from well known Raabe’s criterion for the convergence of series of positive terms, one concludes that ${\sum }_{i=0}^{\infty }{S}_i<\infty$, which means that $\underset{i\to \infty }{lim}{S}_i=S_{dfe}=0$. Note that in Theorem 4, $k_0<1-({\beta }_I{I}_0+{\beta }_E{E}_0)$ is given separately from the general term for the vaccination gain $k_i=\frac{{ϵ}_i}{i^{{\theta }_i}}$ defined for $i>0$ since the gain formula is infinity for $i=0$. □

3. Extended Discrete SEIR Epidemic Model with Delays

An extension of the SEIR model, which includes a finite set of commensurate discrete delays in contrast to the delay-free infection model from Equations (1)–(4), is now discussed. Such a model is the following one:

$$\begin{array}{ccc}\hfill S_{i+1}& =& (1-{\beta }_I{I}_i-{\beta }_E{E}_i)S_i\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill E_{i+1}& =& ({\beta }_I{I}_i+{\beta }_E{E}_i)S_i+E_i-\sum _{j=j_d}^{i-d_1}\left({\mu }_j{E}_j\right)\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill I_{i+1}& =& (1-\gamma )I_i+\sum _{j=j_d}^{i-d_1}\left({\mu }_j{E}_j\right)\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill R_{i+1}& =& \gamma I_i+R_i\hfill \end{array}$$

(25)

for any integer $i\in {\mathbb{Z}}_{0+}$ and any given initial conditions $S_0\ge 0,E_0\ge 0,I_0\ge 0and{R}_0\ge 0$, $j_d=max(i-d-d_1,0)$ with $d\in {\mathbb{Z}}_{+}$, $d_1\in {\mathbb{Z}}_{0+}$. The above model has $d+2$ delays, namely, the one-step-delay of the discretization dynamics from $E_{i+1}$ to $E_i$ plus the contributions of $d+1(\ge 0)$ more extra delays from the sample $i-d_1$ to the sample $i-d-d_1$. The integer $d_1\in {\mathbb{Z}}_{0+}$ reflects a potential extra delay if $d_1>0$, to start to account for extra former delays distinct than the discretization one in the delayed dynamics. If $d_1=0$ then the delayed dynamics starts to influence directly immediately after the one-step discretization step. It becomes technically easier for discussion the mixed normalization and the solution sequence non- negativity to initialize values over $d+1$ consecutive samples, instead for a point initialization for $i=0$ and then to run the model for successive samples. Therefore, the model from Equations (22)–(25) is modified as follows with $d_1=0$, i.e., such that the delayed dynamics starts in the previous sample to the current one and it is affected by a total of d consecutive previous delayed samples with $d\ge 1$:

$$\begin{array}{ccc}\hfill S_{i+1}& =& (1-({\beta }_I{I}_i+{\beta }_E{E}_i))S_i\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill E_{i+1}& =& ({\beta }_I{I}_i+{\beta }_E{E}_i)S_i+E_i-\sum _{j=i-d}^i\left({\mu }_j{E}_j\right)\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill I_{i+1}& =& (1-\gamma )I_i+\sum _{j=i-d}^i\left({\mu }_j{E}_j\right)\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill R_{i+1}& =& \gamma I_i+R_i\hfill \end{array}$$

(29)

for any integer $i(\ge d)\in {\mathbb{Z}}_{+}$ and given $d+1$ sets of initial conditions $S_i\ge 0,E_i\ge 0$, $I_i\ge 0and{R}_i\ge 0$; $i=0,1,\dots ,d$. In the above scheme, the total number of delays is $d+1$ including that of the natural discretization. To distinguish this model from Equations (1)–(4), it is claimed that there is at least one extra delay than the discretization one, i.e., the total minimum number of delays $d+1$ is two since $d\ge 1$. If $d=0$ then the model reduces to Equations (1)–(4) where the only delay is the discretization one. The normalized model from Equations (26)–(29) has a non-negative sequence vector solution provided that all the punctual initial conditions for $i=0,1,\dots ,d$ are normalized as addressed in the next result. Such a result incorporates also some boundedness initial conditions for the infectious and the exposed subpopulations along the relevant number of consecutive delays to fix such initial conditions. Those constraints remember tightly those given in Theorems 1 and 2 for the simpler model from Equations (1)–(4) which possesses just the single one-step discretization delay.

Theorem 5.

Assume that $d_1=0$, ${\beta }_I>0,{\beta }_E>0$, $\gamma \ge {\sum }_{j=0}^d{\mu }_j\ge {\beta }_I+{\beta }_E$, $\gamma \in (0,1]$ and ${\mu }_i\in (0,1]$, $0\le E_i\le \frac{1}{2{\beta }_E}$ and $0\le I_i\le \frac{1}{2{\beta }_I};\forall i\in [1,d]$. Assume also that the initial conditions are non-negative and normalized ($N_0=1$) so that they fulfill the constraints:

$$\left\{\begin{array}{c}min\left(S_i,E_i,I_i,1-(S_i+E_i+I_i)\right)\ge 0\hfill \\ max\left(S_i,E_i,I_i,1-(S_i+E_i+I_i)\right)\le 1\hfill \end{array}\right.;\forall i\in [1,d]$$

(30)

Then the solution is non-negative since:

$$min(S_i,E_i,I_i,R_i)\ge 0,max(S_i,E_i,I_i,R_i)\le 1;i\in {\mathbb{Z}}_{0+}$$

Proof. 

First, note that the constraint $\gamma \ge {\sum }_{j=0}^d{\mu }_j\ge {\beta }_I+{\beta }_E$ implies that $1-{\sum }_{j=0}^d{\mu }_j+{\beta }_I+{\beta }_E\le 1$.

Note also that from Equations (26)–(29) and the non-negativity/normalization constraints of the initial conditions for $j=1$, one has, since ${\sum }_{j=0}^d{\mu }_j\le \gamma \le 1$, that:

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill 0& \hfill {\displaystyle \le \left(1-({\beta }_I+{\beta }_E)\underset{0\le i\le d+j-1}{max}\left(max(I_i,E_i\right)\right)\underset{0\le i\le d+j-1}{min}{S}_i}& \\ \hfill \le S_{d+j}& \hfill {\displaystyle \le \left(1-({\beta }_I+{\beta }_E)\underset{0\le i\le d+j-1}{min}\left(min(I_i,E_i)\right)\right)\underset{0\le i\le d+j-1}{max}{S}_i}& \hfill \le 1\end{array}\end{array}$$

(31)

$$\begin{array}{ccc}\hfill 0& \le & {\displaystyle (1-\sum _{i=j-1}^{j+d-1}{\mu }_i)\underset{i\in [j-1,d+j-1]}{min}{E}_i}\hfill \\ & +& \left({\beta }_I+{\beta }_E\right)\underset{i\in [j-1,d+j-1]}{min}(min(E_i,I_i))\underset{i\in [j-1,d+j-1]}{min}{S}_i\le E_{d+j}\hfill \\ & {\displaystyle \le }& E_{d+j-1}-\left(\sum _{i=j-1}^{j+d-1}{\mu }_i\right)\underset{i\in [j-1,d+j-1]}{min}{E}_i\hfill \\ & +& ({\beta }_E+{\beta }_I)\underset{i\in [j-1,d+j-1]}{max}\left(max(I_i,E_i)\right)\underset{i\in [j-1,d+j-1]}{max}{S}_i\le 1\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill 0& \le I_{d+j}& \le 1-\gamma +\sum _{i=j-1}^{j+d-1}{\mu }_i\le 1\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill 0& \le R_{d+j}& =1-(S_{d+j}+E_{d+j}+i_{d+j})=R_{d+j-1}+\gamma I_{d+j-1}\hfill \end{array}$$

(34)

and by summing up all the subpopulations, one gets that, for $j=1,N_{d+j}=N_{d+1}=N_d=N_0=1$. However, the above constraints also hold via recursive calculations for all $j\in {\mathbb{Z}}_{0+}$ starting from the above equations for $j=0$. The proof of non-negativity is complete. □

Equations (27) and (28) can be jointly described as evolution sequence vectors defined by sets of $d+1$ consecutive samples by a real vector

$${\mathbf{z}}_i={(E_i,E_{i-1},\dots ,E_{i-d},I_i,I_{i-1},\dots ,I_{i-d})}^T$$

through the discrete dynamic system as follows:

$${\mathbf{z}}_{\mathbf{i}+\mathbf{1}}={\mathbf{Az}}_{\mathbf{i}}+{\mathbf{Bu}}_{\mathbf{i}};\forall \mathbf{i}(\ge \mathbf{d})\in {\mathbb{Z}}_{\mathbf{0}+}$$

(35)

subject to a point vector initial condition ${\mathbf{z}}_d\in {\mathbb{R}}^{2(d+1)}$ built from the given d sets of initial conditions from Equations (27) and (28), whose partitioned parametrization is defined as follows:

$$\begin{array}{ccc}& \mathbf{A}=\left[\begin{array}{cc}{\mathbf{A}}_{11}& {\mathbf{0}}_{(d+1)\times (d+1)}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}\end{array}\right]\in {\mathbb{R}}^{2(d+1)\times 2(d+1)};& \mathbf{B}=\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{0}}_{d+1,1}\end{array}\right]\in {\mathbb{R}}^{2(d+1)};\hfill \\ & {\mathbf{u}}_i=\left[\begin{array}{c}({\beta }_I{I}_i+{\beta }_E{E}_i)S_i\\ 0_{2d+1,1}\end{array}\right]\end{array}$$

(36)

with

$$\begin{array}{c}\hfill {\mathbf{A}}_{\mathbf{11}}=\left[\begin{array}{cccc}\mathbf{1}-{\mu }_{\mathbf{i}}& ...& -{\mu }_{\mathbf{i}-\mathbf{d}+\mathbf{1}}& -{\mu }_{\mathbf{i}}-\mathbf{d}\\ & {\mathbf{I}}_{\mathbf{d}}& & {\mathbf{0}}_{\mathbf{d}\times \mathbf{1}}\end{array}\right]\in {\mathbb{R}}^{(\mathbf{d}+\mathbf{1})\times (\mathbf{d}+\mathbf{1})};{\mathbf{B}}_{\mathbf{1}}=\left[\begin{array}{c}\mathbf{1}\\ {\mathbf{0}}_{\mathbf{d}\times \mathbf{1}}\end{array}\right]\in {\mathbb{R}}^{\mathbf{d}+\mathbf{1}}\\ \hfill {\mathbf{A}}_{21}=\left[\begin{array}{ccc}{\mu }_i& ...& {\mu }_{i-d}\\ & {\mathbf{0}}_{d\times (d+1)}\end{array}\right]\in {\mathbb{R}}^{(d+1)\times (d+1)};{\mathbf{A}}_{\mathbf{22}}=\left[\begin{array}{ccc}\mathbf{1}-\gamma \hfill & ...\hfill & {\mathbf{0}}_{\mathbf{1}\times \mathbf{d}}\hfill \\ {\mathbf{I}}_{\mathbf{d}}\hfill & & {\mathbf{0}}_{\mathbf{d}\times \mathbf{1}}\hfill \end{array}\right]\end{array}$$

$I_d$ is the n identity matrix and $0_{n,m}\in {\mathbb{R}}^{n\times m}$ is a zero matrix of order $n\times m$. The recursive solution of Equation (35) becomes:

$${\mathbf{z}}_i={\mathbf{A}}^d{\mathbf{z}}_{i-d}+\sum _{j=i-d}^{i-1}{\mathbf{A}}^{i-j-1}{\mathbf{B}\mathbf{u}}_j$$

(37)

Note that, because of the triangular structure of the matrix A, its $2(d+1)$ eigenvalues are those of $A_{22}$, i.e., $1-\gamma$ plus d zero eigenvalues, and those of $A_{11}$ which, since it is a companion matrix, are the roots of the subsequent discrete polynomial:

$$\eta \left(z\right)=z^{d+1}+({\mu }_i-1)z^d+\sum _{j=i-d}^{i-1}{\mu }_j{z}^{j-i+d}$$

(38)

whose stability property (i.e., its zeros lie in the open unit circle of the complex plane $\left|z\right|<1$) can be investigated through the Jury criterion or via the bilinear transformation $z=\frac{1+s}{1-s}$, which transforms the open unit circle of the complex plane into the open left-hand-side of the complex plane. Thus, the application of the Routh–Hurwitz criterion is useful to investigate if all the zeros of the transformed polynomial $\eta \left(s\right)=\eta (z=\frac{1+s}{1-s})$ belong to the open left-hand-side plane $Re\left(s\right)<0$ or not (see Theorem 2). The next result relies on the convergence of the solution to the disease-free equilibrium point.

Theorem 6.

Assume that $\eta \left(z\right)$ is a strictly stable polynomial and that all the hypotheses of Theorem 5 hold. Then $\underset{i\to \infty }{lim}(1-({\beta }_I{I}_i+{\beta }_E{E}_i)\le 1$, ${\left\{S_i\right\}}_{i=0}^{\infty }\to S_{dfe}\in [0,1]$, ${\left\{R_i\right\}}_{i=0}^{\infty }\to R_{dfe}=(1-S_{dfe})\in [0,1]$ and ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$, ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$. If $\underset{i\to \infty }{lim}(1-{\beta }_I{I}_i-{\beta }_E{E}_i)<1$ then $S_{dfe}=0$ and $liminf{min}_{i\to \infty }(E_i,I_i)\ge 0$.

Proof. 

First note from Theorem 5 that the solution sequence vector is non-negative. Now, note from Equation (26) that ${\left\{S_i\right\}}_{i=0}^{\infty }\to S_{dfe}\in [0,1]$ with $1-{\beta }_I{I}_i-{\beta }_E{E}_i\le 1;\forall i\in {\mathbb{Z}}_{0+}$ and $\underset{i\to \infty }{lim}(1-\beta I_i-{\beta }_E{E}_i)\le 1$, since otherwise, the susceptible subpopulation solution sequence would violate the non-negativity property of Theorem 5. On the other hand, note that:

• If such limit is equal to unity then $S_{dfe}\in [0,1]$ and ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$, ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$ so that ${\left\{R_i\right\}}_{i=0}^{\infty }\to R_{dfe}=1-S_{dfe}$
• If the limit is less than unity then $S_{dfe}=0$ and $liminf{min}_{i\to \infty }(E_i,I_i)\ge 0$. If both ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$, ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$, the proof follows directly. Otherwise, assume that $liminf{min}_{i\to \infty }(E_i,I_i)\ge 0$. Then, since ${\left\{S_i\right\}}_{i=0}^{\infty }\to S_{dfe}=0$ one has from Equation (27) that

  $$\underset{i\to \infty }{lim}\left(E_{i+1}-E_i+\sum _{j=i-d}^i{\mu }_j{E}_j\right)=0$$
  Thus, since $\eta \left(z\right)$ is strictly stable, then ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$. From Equation (28), it follows that ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$ since ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$ and $\gamma <1$. Also ${\left\{R_i\right\}}_{i=0}^{\infty }\to R_{dfe}=1$. Furthermore, ${\left\{S_i\right\}}_{i=0}^{\infty }\to S_{dfe}\in [0,1]$; ${\left\{R_i\right\}}_{i=0}^{\infty }\to R_{dfe}=1-S_{dfe}$ and ${\left\{E_i\right\}}_{i=0}^{\infty }\to 0$, ${\left\{I_i\right\}}_{i=0}^{\infty }\to 0$ so that the solution sequence converges to the disease-free equilibrium point.

□

The subsequent result specifies further Theorem 6 in the sense that it relies on the monotonic convergence to zero levels of the infection in the case that the transmission rates are sufficiently small. In such a way, the levels of the exposed and infectious are strictly decreasing along the transient.

Theorem 7.

Assume that all the hypotheses of Theorem 6 hold. Then $\underset{i\to \infty }{lim}\left|\right|z_i\left|\right|\le \frac{K({\beta }_I+{\beta }_E)}{1-\rho \left(\mathit{A}\right)}$ where $\rho \left(\mathit{A}\right)\in (0,1)$ is the spectral radius of $\mathit{A}$ which is $max(1-\gamma ,{\eta }_M)$, where ${\eta }_M$ is the maximum absolute value of the zeros of $\eta \left(z\right)$ and $K\ge 1$ is some norm-dependent real constant. Furthermore, if $\rho \left(\mathit{A}\right)$ is sufficiently small and ${\beta }_I+{\beta }_E$ is sufficiently small related to $\rho \left(\mathit{A}\right)$ then the sequence ${\displaystyle {\left\{\underset{i-d-1\le j\le i}{max}max(E_j,I_j)\right\}}_{i=0}^{\infty }}$ converges monotonically to zero.

Proof. 

It follows from Equation (26) and the non-negativity of the solution sequences given in Equation (37) that one has, for any vector norm $\parallel u_i\parallel =\left({\beta }_I{I}_i+{\beta }_E{E}_i\right)S_i\le ({\beta }_I+{\beta }_E)\left|\right|z_i\left|\right|$, the following inequality:

$$\begin{array}{ccc}\hfill \parallel z_i\parallel & \le & \parallel {\mathbf{A}}_d\parallel \parallel z_{i-d}\parallel +\sum _{j=i-d}^{i-1}\parallel {\mathbf{A}}^{i-j-1}\parallel \parallel u_j\parallel \hfill \\ & \le & \parallel {\mathbf{A}}_d\parallel \parallel z_{i-d}\parallel +\left(\sum _{j=i-d}^{i-1}\parallel {\mathbf{A}}^{i-j-1}\parallel \right)\underset{j\in [i-d,i-1]}{sup}\parallel u_j\parallel \hfill \\ & \le & \left[\parallel {\mathbf{A}}_d\parallel +\underset{j\in [i-d,i-1]}{sup}\parallel z_j\parallel ({\beta }_I+{\beta }_E)\left(\sum _{j=d}^{i-1}\parallel {\mathbf{A}}^{i-j-1}\parallel \right)\right]\underset{j\in [i-d,i-1]}{sup}\parallel z_j\parallel \hfill \\ & \le & \left[K\rho \left({\mathbf{A}}_d\right)+\underset{j\in [i-d,i-1]}{sup}\parallel z_j\parallel \frac{K\left(1-{\rho }^{i-d}\left(\mathbf{A}\right)\right)}{1-\rho \left(\mathbf{A}\right)}({\beta }_I+{\beta }_E)\right]\underset{j\in [i-d,i-1]}{sup}\parallel z_j\parallel \hfill \\ & \le & K\rho {\left(\mathbf{A}\right)}^{i-d}\underset{j\in [i-d,i-1]}{sup}\parallel z_j\parallel +\frac{K}{1-\rho \left(\mathbf{A}\right)}({\beta }_I+{\beta }_E)\underset{j\in [i-d,i-1]}{sup}{\parallel z_j\parallel }^2\hfill \end{array}$$

(39)

since $\parallel \mathbf{A}\parallel \le K\rho \left(\mathbf{A}\right)$ and the fact that the model is normalized implies from the definition of $u_i$ in Equation (36) that:

$$\frac{K}{1-\rho \left(\mathbf{A}\right)}({\beta }_I+{\beta }_E)\underset{d\le j\le i-1}{sup}{\parallel z_j\parallel }^2\le \frac{K}{1-\rho \left(\mathbf{A}\right)}({\beta }_I+{\beta }_E)$$

(40)

Now, one gets from Equation (39) that $\underset{i\to \infty }{lim}\parallel z_i\parallel \le \frac{K({\beta }_I+{\beta }_E)}{1-\rho \left(\mathbf{A}\right)}$ since $\underset{i\to \infty }{lim}\rho {\left(\mathbf{A}\right)}^{i-d}=0$. Furthermore, one has from Equations (35) and (36) that

$$\parallel z_{i+1}\parallel \le \parallel \mathbf{A}\parallel \parallel z_i\parallel +\parallel \mathbf{B}\parallel \parallel {\mathbf{u}}_i\parallel \le K\rho \left(\mathbf{A}\right)\parallel z_i\parallel +({\beta }_I+{\beta }_E)S_{i}max(E_i,I_i)$$

(41)

Assume that the infinity vector and vector-matrix induced norms are used in the above formula so that

$$\underset{i-d\le j\le i}{max}max\left(E_{j+1},I_{j+1}\right)\le \left(K\rho \left(\mathbf{A}\right)+({\beta }_I+{\beta }_E)\right)\underset{i-d\le j\le i}{max}max(E_j,I_j)$$

(42)

if ${\beta }_I+{\beta }_E<1-K\rho \left(\mathbf{A}\right)$, then the sequence $\left\{{max}_{i-d\le j\le i}max(E_j,I_j)\right\}$ converges monotonically to zero. □

It is of interest to have some testing tool to exclude eventual non-trivial oscillations of the solutions of the above model under certain model parametrizations [53]. One sees from Equations (22)–(29) that ${\left\{S_i\right\}}_{i=0}^{\infty }$ converges asymptotically to a limit $S^{*}\ge 0$ and ${\left\{R_i\right\}}_{i=0}^{\infty }$ can only converge to a limit oscillation if ${\left\{I_i\right\}}_{i=0}^{\infty }$ converges asymptotically to an oscillation. Therefore, the whole system (22)–(25) converges to a limit cycle if and only if ${\left\{E_i\right\}}_{i=0}^{\infty }$ and ${\left\{I_i\right\}}_{i=0}^{\infty }$ converge to limit oscillations. Thus for $S\equiv S^{*}$ the limit subsystem (22)–(25) can be described through the following auxiliary discrete system got “ad hoc” for this case from (27)–(28):

$$z_{i+1}=\mathbf{A}{z}_i+S^{*}({\beta }_I{I}_i+{\beta }_E{E}_i){\mathbf{e}}_1;\forall i(\ge d)\in {\mathbb{Z}}_{0+}$$

(43)

where ${\mathbf{e}}_1\in {\mathbb{R}}^{2(d+1)}$ is the unity Euclidean vector of first component unity and, if there is a non-trivial limit oscillation solution of Equation (43) then ${\displaystyle \underset{i\to \infty }{lim}{z}_i=z^{*}=\left[\begin{array}{c}{\lambda }^{*}\\ {\omega }^{*}\end{array}\right]\ne 0\left(\in {\mathbb{R}}^{2(d+1)}\right)}$ which satisfies from (43) the following constraint:

$$\left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-S^{*}{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]{\mathbf{z}}^{*}=0$$

(44)

so that $z^{*}\in Ker\left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-S^{*}{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]$ and then there is no non-trivial limit cycle if

$$\left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-S^{*}{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]$$

is non-singular since then $z^{*}=0\left(\in {\mathbb{R}}^{2(d+1)}\right)$. Since ${\left\{S_i\right\}}_{i=0}^{\infty }$ is non-increasing then $S^{*}=S_0-{\tilde{S}}^{*}$ for some ${\tilde{S}}^{*}={\tilde{S}}^{*}\left(S_0\right)$. Then, from Banach’s Perturbation Lemma [54], one gets

$$\begin{array}{c}\hfill \left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-S^{*}{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]\\ \hfill =\left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-(S_0-{\tilde{S}}^{*}){\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]\end{array}$$

(45)

is non-singular if $\left[{\mathbf{I}}_{2(d+1)}-\mathbf{A}-S_0{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]$ is non-singular and ${\tilde{S}}^{*}({\beta }_I+{\beta }_E)<{∥{\left[I_{2(d+1)}-\mathbf{A}-S_0{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T\right]}^{-1}∥}_1$ since then $\parallel {\tilde{S}}^{*}{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T){\parallel }_1\le {\tilde{S}}^{*}({\beta }_I+{\beta }_E)<{∥{\left[I_{2(d+1)}-\mathbf{A}-S_0{\mathbf{e}}_1({\beta }_I{\mathbf{e}}_1^T+{\beta }_E{\mathbf{e}}_{d+2}^T)\right]}^{-1}∥}_1$.

4. Data Processing

We will now proceed to determine the values of the transmission rates ${\beta }_I$ and ${\beta }_E$ for the SEIR models using the available data in the cited countries.

4.1. Classification of Control Strategies

Except for Sweden [41], the countries presented in this paper (i.e., Spain [34], Italy [35], France [36], Portugal [39], United Kingdom [37], Germany [38] and Norway [40]) have implemented different strategies in order to change the course of the epidemic at all levels. There are many different ambient conditions that may affect the value of both infectivity rates, such as the weather conditions or social awareness to the disease, but we will assume that the greater the level of lockdown that the country submits itself to will mean more influence over the transmission rates. Thus, with the purpose of having the same criteria independently of the politics of each country, the following classification of the diverse isolation measures applied to the population is made:

• Stage 1: Total lockdown (Blue background in the graphics) People must stay at home unless necessary. Going out for vital activities such as work, groceries and exercise may be allowed.
• Stage 2: Outdoor lockdown (Red background in the graphics) Necessary tasks and outdoor contact allowed. Going out is permitted with established schedules. In addition, schools are open, family visits allowed, etc. Activities in closed spaces with moderate risk of transmission (bars, restaurants, museums, gyms, cinemas, etc.) are shut down.
• Stage 3: Indoor lockdown (Yellow background in the graphics) Public indoor events allowed. Activities in closed spaces are opened with certain restrictions, such as a reduction in the schedules and/or the maximum capacity of the locations.
• Stage 4: New normality (White background in the graphics) All services are available although they may require basic hygiene and prophylactic measures. Borders are re-opened with routine health check-ups applied to migrants.

As we can see, the social distancing measurements adopted at each phase present a relaxation of those of the previous ones gradually loosening the isolation up to a normal situation.

4.2. Processing the Available Data

In order to calculate the transmission rates we will make some assumptions:

• During the lockdown stages adopted in the different studied countries, the social interaction of an exposed individual is similar to a diagnosed infectious one.
• The exposed transition rate during the lockdown stage will be equal to the infectious transition rate ${\beta }_{E_0}={\beta }_{I_0}$.
• The diagnosis methods and the prophylactic measures around these diagnosed infectious individuals will improve after the beginning of the epidemic.
• The value of ${\beta }_I$ will be reduced partially as the awareness of the disease and the risk around infected individuals are reduced because of the social distancing strategies adopted in the country.
• The value of the transmission rate of the diagnosed infectious ${\beta }_I$ after the initial lockdown will remain constant as the official approach to a diagnosed individual will remain the same independently of the lockdown stage.
• The bias in the influence of the non-measured infected individuals will be constant through all the time.

The ${\beta }_I$ and ${\beta }_E$ will be calculated then in two steps. First, we will take the available data from the lockdown stages applied on a population, considering equal the initial transmission rates ${\beta }_{E_0}={\beta }_{I_0}$. The value of the total cumulative cases, which will include all the infected subpopulations, is the same independently of the SEIR model chosen from Equations (1)–(4) or (26)–(29), namely, it is $C{C}_i=E_i+I_i+R_i$. We obtain that

$$\Delta C{C}_i=\Delta E_i+\Delta I_i+\Delta R_i=S_i({\beta }_I{I}_i+{\beta }_E{E}_i)={\beta }_{I_0}{S}_i(I_i+E_i)$$

(46)

Since ${\beta }_E={\beta }_I$ in this first step Then, as the value for the susceptible at any day i would be $S_i=1-C{C}_i$, we will get the equation for ${\beta }_I$ and ${\beta }_E$ for the lockdown time

$${\beta }_{E_0}={\beta }_{I_0}=\frac{\Delta C{C}_i}{(I_i+E_i)(1-C{C}_i)}$$

(47)

A linear regression will be used to obtain the ${\beta }_{I_0}$ at the initial total lockdown stage. Then, the transmission rate for the infectious at the next stages of the lockdown ${\beta }_I$ will be reduced to a fraction of the original value as the prophylactic measures and the awareness of the disease increases. Then, the ${\beta }_I$ for the rest of the lockdown stages will be ${\beta }_I=k{\beta }_{I_0}$ with the factor $k\in [0,1]$ set to a point as the transmission rate of the exposed ${\beta }_E$ is calculated. The ${\beta }_E$ will be obtained applying this k to the rest of the stages and calculating again through linear regression using the data of the days of each stage as

$${\beta }_E=\frac{\Delta C{C}_i}{E_i(1-C{C}_i)}-\frac{{\beta }_I{I}_i}{E_I}$$

(48)

Different k will be studied in order to fit the most appropriate ${\beta }_E$.

4.2.1. Estimation of the Exposed Subpopulation

As we can see in Equations (46)–(48) and Equations (1)–(4) and (26)–(29) respectively, when calculating the values of ${\beta }_I$ and ${\beta }_E$, the value of the exposed subpopulation is as important as the infectious one. However, using only the official repositories, in which the different authorities share their data related to COVID-19, we cannot obtain the number of individuals incubating the virus at any given time, as it is obviously something which is not directly observed and thus not registered. In order to extrapolate such number from the available data for the extended model, we will assume:

• That each infected individual comes from a exposed individual.
• That each exposed individual will eventually be registered as infected. This assumes that they will be tested due to developing symptoms or being in close contact to someone who has.
• That each exposed individual will become infected individual in a range of days from $k_0$ to $k_1$, with a probability of transition based on a normal distribution around a central value [55].
• The density function of the probability of transition will be defined from experimental observations [56,57] by the mean value (4 days) an the variance (3 days), as seen as in [58].

The estimation of the exposed subpopulation will be done for all dates in a single step, using the known values of new cases of all dates. Now, let $E_k$ and $C_k$ be the number of exposed individuals and the new infected cases at the day k, respectively. Then, the new infected individuals at day k will depend on the probability of transition and the number of exposed individuals on previous days. In this case, as:

$$C_k=g_{k_0}{E}_{k-k_0}+\dots +g_{k_1}{E}_{k-k_1}=\sum _{m=k_0}^{k_1}{g}_m{E}_{k-m}$$

(49)

with the fraction of individuals that transient to the infected individuals $g_i=\frac{{\phi }_{\mu ,\sigma }\left(i\right)}{{\sum }_{j=k_0}^{j=k_1}{\phi }_{\mu ,\sigma }\left(j\right)}$, being ${\phi }_{\mu ,\sigma }\left(x\right)$ the value of the probability density function of a normal distribution with mean $\mu$ and standard deviation $\sigma$ at specific values $i\in [k_0,k_0+1,\dots ,k_1]$. Observe that these values of $g_k$ correspond to a special case of the extended model in which $\sum {\mu }_i=1$. Then we will build a database with n data and a $(n-k_0)\times (n-k_0)$ matrix A with the values $g_{k_0},\dots ,g_{k_1}$ that multiplying the vector of historical exposed $\mathbf{E}={(E_1,\dots ,E_{n-k_0})}^T$ we would obtain the vector of new cases with a delay of $k_0$ days, $\mathbf{C}={(C_{k_0+1},\dots ,C_n)}^T$. As the known vector is $\mathbf{C}$ we will invert the matrix A, so that multiplying its inverse to the right to the vector $\mathbf{C}$, we get the desired vector $\mathbf{E}$ as described in the following equation

$$\left(\begin{array}{c}{C}_{k_0+1}\\ ⋮\\ C_{n-1}\\ C_n\end{array}\right)=\mathbf{A}\left(\begin{array}{c}{E}_1\\ E_2\\ ⋮\\ E_{n-k_0}\end{array}\right)⟹\left(\begin{array}{c}{E}_1\\ E_2\\ ⋮\\ E_{n-k_0}\end{array}\right)={\mathbf{A}}^{-1}\left(\begin{array}{c}{C}_{k_0+1}\\ ⋮\\ C_{n-1}\\ C_n\end{array}\right)$$

(50)

In this approximation, we will assume that there is no exposed individuals before day 1, so the infected individuals for the first days must come from these initial transitions. This assumption is reasonable because the number of cases at these moments is relatively low, so the value of exposed individuals must also be low. Consequently, the transitions of the first rows from 1 to $k_1-k_0$ will be normalized, so that every row in A equals to 1. This way, we define the matrix $\mathbf{A}=<a_{i,j}>$ with:

$$\begin{array}{c}\hfill a_{i,j}=\left\{\begin{array}{c}i\le (k_1-k_0)+1:\left\{\begin{array}{cc}j\le i\hfill & :\frac{g_{k_0+i-j}}{{\sum }_{l_0}^i{g}_{k_0+l}}\hfill \\ j>i\hfill & :0\hfill \end{array}\right.\hfill \\ i>(k_1-k_0)+1:\left\{\begin{array}{cc}j\in (0,i+k_0-k_1)\cup (i,n-k_0]\hfill & :0\hfill \\ j\in [i,i+k_0-k_1]\hfill & :g_{k_0+i-j}\hfill \end{array}\right.\hfill \end{array}\right.\end{array}$$

(51)

As matrix A is lower triangular, ${\mathbf{A}}^{-1}$ will be lower triangular too. Finally, the curve of the values of the exposed individuals over time will be smoothed through weighted moving average [59], in order to filter out the possible noises in measurements and bureaucratic errors when publishing the results. We can see that for the model from (1)–(4), that the extrapolation of the exposed population is trivial, as it can be easily extracted from $E_i=C_i/\mu$, being $C_i$ the daily value of the new cases.

4.2.2. Estimation of the Recovered Subpopulation

The recovered subpopulation (R), however, will be obtained from the average ratio of recovery from the disease [56,57], as shown in the transition of Equations (3) and (4). Observe that, due to share a similar dynamics in their non-interaction with the susceptible or infected population, the recovered individuals and the individuals that die from the disease are both considered as recovered.

4.2.3. Summer Phase

Observe that the process of infection is an stochastic process, which is specially relevant during the summer period, as during this season the percentage of daily cases of new infected over the total population is quite small. The calculation about such a variability in the value of new infected individuals over such small time cannot be used in a useful way. The transmission rate obtained during the initial stage and during the summer period shown in the figures will not present the required minimum of new infections per number of individuals that belong to the population. Thus the rates ${\beta }_I$ and ${\beta }_E$ will vary too much in a very small period of time, so the average value will not be reliable. These stages will be shown at the figures in a grey background in contrast to the other stages.

5. Numerical Results of the SEIR Model

In this section we will calculate the parameters ${\beta }_I$ and ${\beta }_E$ of different countries in Europe via multiple linear regression, adjusted to the SEIR model from Equations (1)–(4) and (26)–(29) respectively. Data from the official health-care institutions from countries of Spain, Italy, United Kingdom, Germany, Portugal, France, Norway and Sweden will be used. The figures are completed with the new daily cases per 100,000 inhabitants in such a way the reader can observe its evolution. In every case, the curves describing the cumulative and direct values of the infected and exposed individuals are smoothed through weighted moving average [59], in order to filter out the bureaucratic errors when publishing the results and other possible noises in measurement. The analysis period will be from the beginning of March 2020 to the end of February 2021. Each background color of the graphs represents a different stage of the social distancing measures adopted by the territorial government.

5.1. Procedure to Calculate ${\beta }_E$ and ${\beta }_I$

First the different NPI measures of control of COVID-19 at different time intervals will be set, and the average ${\beta }_I$ and ${\beta }_E$ of the different countries, with the two SEIR models will be calculated as described in Section 4.2. The average values will be calculated globally for each country, and then at a lower level of territorial organization (i.e., provinces, regions or counties), when the data is available.

5.2. Determination of the k Factor

We will first determinate the factor k from the equation for the reduction of the infectious transmission constant after the initial lockdown ${\beta }_I=k{\beta }_{I_0}$. This factor will take into account the general change of attitude when acting towards a person diagnosed with COVID-19, as the initial prophylactic measures and prevention methods, such as the face masks, were not sufficiently popularized on the society compared to the actual implementations, which have their positive impact [60,61,62]. We have chosen first France in order to calibrate the best k that fits the results for the expected value of ${\beta }_E$, and then we use this value of k in the rest of the countries.

5.3. France

Data from the official health-care institutions of France is used to calculate the average transition rates during the first lockdown. In this case, the available data did not contain the COVID-19 cases distributed regionally, so only the figures showing the data in the whole country are displayed. In

Table 1. Data for France: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
17 March 2020	7730	11.2	37.6	225
11 May 2020	139,519	16.7	0.5	26,643
1 August 2020	189,547	20.7	1.9	30,251
11 September 2020	363,350	143.5	1.9	30,893
30 October 2020	1,331,984	741.4	2.2	36,565
28 November 2020	2,208,699	378.9	0.5	52,127
12 February 2021	3,427,386	408.4	1.0	81,448

it can be seen the numbers of cases in France at the days the different intervals start and end, as well as different indicators of infection and number of deaths related to COVID-19. The total CFR (Case fatality rate) for France, i.e., the percentage of death individuals per diagnosis is 2.58%.

5.3.1. Extended SEIR Model

We can see in Figure 2 the average ${\beta }_E$ calculated through the lockdown stages in France. The transmission rate related to the individuals diagnosed as Infected, ${\beta }_{I_0}$ will be, as said in Section 4.2 the same as the ${\beta }_{E_0}$ during the lockdown stage. The exposed subpopulation will be obtained from the extrapolation of the newly infected as described in Equations (27), (28), (49) and (50). Then, due to the popularization of diverse NPI, the infectious transmission rate will be reduced by a factor k, for the relation ${\beta }_I=k{\beta }_{I_0}$, which we will estimate through calibration. A sweep of k between values 0.5 to 1 is made and the average ${\beta }_E$ calculated through lineal regression from Equation (48).

5.3.2. $\mu$-SEIR Model

Again, we will extrapolate the exposed subpopulation from the model Equations (1)–(4), and recalculate ${\beta }_{E_0}$ and ${\beta }_{I_0}$ of the first lockdown Stage. Then, as in the previous figure, ${\beta }_E$ will be calculated given different stages.

We estimate that the different possible values of ${\beta }_E$, from Figure 2 and Figure 3, obtained by using a k-factor $k=0.7$, as it presents a behavior more adequate to the expected value of ${\beta }_E$. Once this value is established, it will be used in the rest of the countries in Europe, as the improvement on the treatment of infected individuals once the first wave passed is assumed the same.

5.4. Norway

Information about the impact of the disease in Norway at the points of transition from the different stages can be seen in

Table 2. Data for Norway: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
12 March 2020	860	15.9	15.9	0
20 April 2020	7106	25.1	0.4	154
8 September 2020	11,448	20.9	1.8	264
5 November 2020	22,467	104.2	3.2	284
12 February 2021	65,457	65.6	0.9	592

. The total case fatality rate of Norway is 0.9%.

5.4.1. Extended SEIR Model

As in the SEIR model with normal distribution for France, we calculate the exposed subpopulation as said in Section 4.2 and calculate the transmission rates ${\beta }_I$ and ${\beta }_E$, shown in Figure 4.

5.4.2. $\mu$-SEIR Model

Again, the transmission rates ${\beta }_I$ and ${\beta }_E$ are calculated for the SEIR model from Equations (1)–(4) and shown in Figure 5.

5.5. Sweden

Sweden is a special case in Europe, as they did not consider until the beginning of November that any special lockdown were needed, so the first wave of COVID-19 presented an incidence much higher than their neighboring country Norway, as we see in Table 2 and

Table 3. Data for Sweden: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
23 March 2020	2169	17.8	5.7	33
21 July 2020	74,766	35.1	0.4	5670
5 September 2020	85,500	24.4	0.7	5848
10 November 2020	166,956	470.7	1.6	6092
12 February 2021	609,306	400.1	1.0	12,449

comparing the incidence rates (green dotted line) and the infectious rates. The data regarding CFR of Sweden was 2.19%.

5.5.1. Extended SEIR Model

The ${\beta }_I$ from the Norway normal distribution model is used due to the similar social structures and interaction frequency [63]. We proceed as in previous countries in order to calculate the ${\beta }_E$ during the first wave and the late social distancing measures. Figure 6 displays the results.

5.5.2. $\mu$-SEIR Model

The ${\beta }_I$ from the Norway in the $\mu$-SEIR model is used as in the previous section due to the same reasons, and the ${\beta }_E$ is calculated again for the first wave of the disease and the late social distancing measures. The results are displayed in Figure 7

5.6. Portugal

In the case of Portugal, the data of incidence on the affected population at the starting point of the lockdown stages is presented in

Table 4. Data for Portugal: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
23 March 2020	642	6.2	106	1
4 May 2020	25,524	45.3	0.6	1063
1 July 2020	42,523	47.1	1.2	1579
15 September 2020	65,021	65.8	1.9	1875
15 January 2021	528,469	1047.4	2.2	8543
12 February 2021	781,223	802.6	0.5	15,034

. The total case fatality rate of Portugal up to February 2021 is 2.5%.

5.6.1. Extended SEIR Model

As in previous countries the ${\beta }_E$ is calculated for the different stages of the disease and shown in Figure 8. Observe that after the “new normality” the transmission rate for the exposed rises unexpectedly.

5.6.2. $\mu$-SEIR Model

The transmission rate for the exposed population in this model is similar to the previous one as seen in Figure 9, in that the second lockdown measures do not seem to affect as much as the first one.

5.7. United Kingdom

The United Kingdom had the second worst case fatality rate of all the countries here shown, with a value of 3.2%. The strategy adopted by the country up to February 2021 is to implement full lockdown measures The incidence data at the beginning of the different lockdown stages are presented in

Table 5. Data for the UK: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of change of I in 14 days, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
23 March 2020	12,639	18	19.9	938
1 June2020	257,702	37.9	0.7	38,294
20 August 2020	325,468	22.37	1.4	41,494
5 November 2020	1,183,510	464.9	1.3	49,271
2 December 2020	1,696,375	317.5	0.7	61,050
4 January 2021	2,880,767	1042.7	1.8	79,323
12 February 2021	4,040,037	328.2	0.6	118,097

.

5.7.1. Extended SEIR Model

We calculate ${\beta }_E$ as in the previous countries. Observe in Figure 10 the drop of the values at the closed lockdown stages (Blue background zone).

In this case, we can also use the linear regression technique for the individual 64 different counties that constitute the UK. A series of histograms related to each stage can be seen at Figure 11 and their associated statistical data at

Table 6. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in the UK.

Phase	Mean	Variance	CI 95%	p
Lockdown	0.08	$10\times {10}^{-4}$	(80.7, 82.5) $\times {10}^{-3}$
Indoor contact	0.46	12 $\times {10}^{-4}$	(452.6, 469.5) $\times {10}^{-3}$	0
Second lockdown	0.22	134 $\times {10}^{-4}$	(192.0, 248.6) $\times {10}^{-3}$	0
Indoor contact	0.47	80 $\times {10}^{-4}$	(452.0, 495.7) $\times {10}^{-3}$	0
Third lockdown	0.17	69 $\times {10}^{-4}$	(144.8, 185.4) $\times {10}^{-3}$	0

. There, the mean value of the ${\beta }_E$s, their variance and the confident interval are displayed, as well as the p-value between histograms to show the notable impact of the strategies on the different stages.

5.7.2. $\mu$-SEIR Model

Again, the transmission rate for the exposed subpopulation ${\beta }_E$ is calculated and shown in Figure 12.

The histograms of the different stages describing the ${\beta }_E$s of the different counties of the UK are presented in Figure 13, slightly different to the previous one.

In

Table 7. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in the UK.

Phase	Mean	Variance	CI 95%	p
Lockdown	$64.8\times {10}^{-3}$	$10\times {10}^{-4}$	(64, 66) $\times {10}^{-3}$
Indoor contact	$169\times {10}^{-3}$	$3.2\times {10}^{-4}$	(164, 173) $\times {10}^{-3}$	0
Second lockdown	$83\times {10}^{-3}$	$14.3\times {10}^{-4}$	(74, 93) $\times {10}^{-3}$	0
Indoor contact	$186\times {10}^{-3}$	$10.2\times {10}^{-4}$	(178, 194) $\times {10}^{-3}$	0
Third lockdown	$5.3\times {10}^{-3}$	$5.4\times {10}^{-4}$	(48, 59) $\times {10}^{-3}$	0

, the statistical data associated to the histograms is seen.

5.8. Germany

For Germany, the case fatality rate is 3.1%. The incidence of COVID-19 at the different starting points of the stages of lockdown can be seen in

Table 8. Data for Germany: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
23 March 2020	35,985	40.8	16.3	708
20 April 2020	146,504	50.5	0.6	7537
14 September 2020	264,405	23.1	1.2	9593
2 November 2020	579,340	234.8	2.7	13,417
24 December 2020	1,614,002	397.52	1.3	43,701
10 January 2021	1,933,486	315.9	0.9	55,945
12 February 2021	2,326,814	147.8	0.7	69,327

.

5.8.1. Extended SEIR Model

A calculation of the ${\beta }_E$ is made at the different lockdown stages. Observe in Figure 14 a diminishing of the value of ${\beta }_E$ in the Stage 2 after the second Stage 1 lockdown, whose value is almost identical to the initial value of ${\beta }_E$ at the first Stage 1 of lockdown in the spring of 2020.

The impact of the lockdown stages on the 33 different regions of Germany can be easily seen at the histograms of Figure 15.

The data and the significant statistical differences between the stages shown with the p-test are displayed at

Table 9. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Germany.

Phase	Mean	Variance	CI 95%	p
Lockdown	$7.9\times {10}^{-2}$	$3.0\times {10}^{-3}$	(78, 79) $\times {10}^{-3}$
Indoor contact	$56.8\times {10}^{-2}$	$1.0\times {10}^{-3}$	(557, 579) $\times {10}^{-3}$	0
Outdoor contact	$38.1\times {10}^{-2}$	$5.8\times {10}^{-3}$	(355, 407) $\times {10}^{-3}$	0
Second Lockdown	$21.4\times {10}^{-2}$	$12.5\times {10}^{-3}$	(176, 252) $\times {10}^{-3}$	<${10}^{-6}$
Outdoor contact	$8.2\times {10}^{-2}$	$7.3\times {10}^{-3}$	(53, 111) $\times {10}^{-3}$	<${10}^{-6}$

.

5.8.2. $\mu$-SEIR Model

A second calculation of the transmission rate ${\beta }_E$ is made for the SEIR model from Equations (1)–(4). As in the previous model we can see in Figure 16 how the ${\beta }_E$ decreases as the lockdown intensifies, except for the final Stage.

The histograms shown at Figure 17 and their related data shown in

Table 10. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Germany.

Phase	Mean	Variance	CI 95%	p
Lockdown	$6.0\times {10}^{-2}$	$0.4\times {10}^{-4}$	(5.9, 6.0) $\times {10}^{-2}$
Indoor contact	$22.2\times {10}^{-2}$	$1.6\times {10}^{-4}$	(21.7, 22.6) $\times {10}^{-2}$	0
Outdoor contact	$13.7\times {10}^{-2}$	$6.1\times {10}^{-4}$	(12.8, 14.5) $\times {10}^{-2}$	0
Second Lockdown	$9.9\times {10}^{-2}$	$10.1\times {10}^{-4}$	(8.8, 11.0) $\times {10}^{-2}$	<${10}^{-8}$
Outdoor contact	$5.2\times {10}^{-2}$	$6.6\times {10}^{-4}$	(4.3, 6.1) $\times {10}^{-2}$	<${10}^{-8}$

prove a clear impact and a significant statistical difference on the ${\beta }_E$ between the lockdown stages.

5.9. Spain

For Spain, the data shows that the case fatality rate is 2.3%. The incidence of COVID-19 at the different starting points of the stages of lockdown can be seen in

Table 11. Data for Spain: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
15 March 2020	12,874	26.9	108.9	508
4 May 2020	228,788	44.9	0.5	26,878
10 July 2020	266,181	16.6	1.8	29,796
15 August 2020	383,287	145.6	1.9	30,245
6 November 2020	1,453,433	619.4	1.4	41,167
4 December 2020	1,781,385	246.2	0.6	48,580
12 February 2021	3,106,326	471.2	0.5	66,858

.

5.9.1. Extended SEIR Model

A calculation of the ${\beta }_E$ is made at the different lockdown stages as described in Section 4.2. The values of the transmission parameter ${\beta }_E$ are shown at Figure 18.

A series of histograms related to the transmission rates of the exposed individually calculated at the 50 provinces of Spain at the different stages of the lockdown is shown at Figure 19.

The data regarding the histograms is displayed in

Table 12. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Spain.

Phase	Mean	Variance	CI 95%	p
Lockdown	$7.8\times {10}^{-2}$	$0.2\times {10}^{-4}$	(7.6, 7.9) $\times {10}^{-2}$
New normality	$57.1\times {10}^{-2}$	$37.1\times {10}^{-4}$	(55.5, 58.9) $\times {10}^{-2}$	0
Indoor contact	$43.3\times {10}^{-2}$	$29.1\times {10}^{-4}$	(41.8, 44.8) $\times {10}^{-2}$	0
Outdoor contact	$15.7\times {10}^{-2}$	$201.1\times {10}^{-4}$	(11.7, 19.6) $\times {10}^{-2}$	0
Indoor contact	$45.4\times {10}^{-2}$	$174.7\times {10}^{-4}$	(41.7, 49.0) $\times {10}^{-2}$	0

, showing the mean value of the ${\beta }_E$s, their variance, the confident interval and the p-test done between the stages of lockdown to check their statistical differences.

5.9.2. $\mu$-SEIR Model

As in the previous SEIR model, the transmission rates for the exposed is calculated at the different lockdown stages as shown in Figure 20.

Again, a series of histograms is made for the ${\beta }_E$s of the 50 provinces of Spain (we have taken out city provinces Ceuta and Melilla) to check for the impact of the lockdown and social distancing measures on the transmission rates as seen in Figure 21.

The data regarding their mean values, their variance, confidence intervals and the p-test between the histograms can be seen at

Table 13. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Spain.

Phase	Mean	Variance	CI 95%	p
Lockdown	$6.1\times {10}^{-2}$	$0.1\times {10}^{-4}$	(6.0, 6.1) $\times {10}^{-2}$
New normality	$24.5\times {10}^{-2}$	$27.6\times {10}^{-4}$	(23.0, 25.9) $\times {10}^{-2}$	0
Indoor contact	$15.4\times {10}^{-2}$	$4.4\times {10}^{-4}$	(14.9, 16.0) $\times {10}^{-2}$	0
Outdoor contact	$5.1\times {10}^{-2}$	$16.6\times {10}^{-4}$	(4.0, 6.2) $\times {10}^{-2}$	0
Indoor contact	$17.3\times {10}^{-2}$	$39.2\times {10}^{-4}$	(15.6, 19.0) $\times {10}^{-2}$	0

.

5.10. Italy

Italy shows the greatest case fatality rate of all the countries studied in this paper: 3.7 deaths per 100 cases on average. The incidence of COVID-19 at the different starting points of the stages of lockdown can be seen in

Table 14. Data for Italy: Total cases, Cases per 100,000 inhabitants in the last 14 days, rate of Cases in the last 14 days and the previous 14 days interval, Death related to COVID-19.

Date	Total Cases	I	$\frac{{\mathit{I}}_{\mathit{i}}}{{\mathit{I}}_{\mathit{i}-14}}$	Death
10 March 2020	10,149	16.3	31.1	631
4 May 2020	211,938	50.9	0.7	29,079
27 August 2020	263,949	19.4	2.5	35,463
25 October 2020	525,782	283.2	3.9	37,338
24 December 2020	2,009,317	368.3	0.9	70,900
6 January 2021	2,201,945	349.3	1.0	76,877
12 February 2021	2,697,296	278.9	1.1	93,045

.

5.10.1. Extended SEIR Model

The calculation of the ${\beta }_E$ from the SEIR model described in Equations (26)–(29) is made as in previous countries and the results are shown in Figure 22.

Again, a series of histograms is made in order to describe the distribution of ${\beta }_E$s calculated for each of the 107 regions of Italy at the different lockdown stages which have been applied since early 2020 as seen in Figure 23.

The statistical information acquired from these histograms is displayed at

Table 15. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Italy.

Phase	Mean	Variance	CI 95%	p
Lockdown	0.079	$0.3\times {10}^{-4}$	(7.8, 8.0) $\times {10}^{-2}$
New normality	0.568	$21.5\times {10}^{-4}$	(55.9, 57.6) $\times {10}^{-2}$	0
Indoor contact	0.399	$24.5\times {10}^{-4}$	(39.0, 40.9) $\times {10}^{-2}$	0
Second lockdown	0.206	$226.2\times {10}^{-4}$	(17.7, 23.4) $\times {10}^{-2}$	0
Indoor contact	0.268	$154.1\times {10}^{-4}$	(24.5, 29.2) $\times {10}^{-2}$	<${10}^{-3}$

.

5.10.2. $\mu$-SEIR Model

Calculations of the transmission rates of the exposed subpopulation in the different lockdown stages is made in Figure 24.

A series of histogram is made as in the previous cases for the 107 regions of Italy at the different lockdown stages, as seen in Figure 25.

Furthermore, the data regarding these histograms is displayed in

Table 16. Histogram for the mean values of ${\beta }_E$ at the different stages of lockdown in Italy.

Phase	Mean	Variance	CI 95%	p
Lockdown	$6.2\times {10}^{-2}$	$0.2\times {10}^{-4}$	(6.1, 6.2) $\times {10}^{-2}$
New normality	$26.3\times {10}^{-2}$	$10.3\times {10}^{-4}$	(25.6, 26.9) $\times {10}^{-2}$	0
Indoor contact	$13.0\times {10}^{-2}$	$1.7\times {10}^{-4}$	(12.7, 13.2) $\times {10}^{-2}$	0
Second lockdown	$11.4\times {10}^{-2}$	$27.0\times {10}^{-4}$	(10.4, 12.4) $\times {10}^{-2}$	$3\times {10}^{-3}$
Indoor contact	$10.2\times {10}^{-2}$	$14.2\times {10}^{-4}$	(9.5, 11.0) $\times {10}^{-2}$	$6\times {10}^{-4}$

.

6. Discussion

6.1. France and Portugal

As said in the previous section, even though France and Portugal are big enough to acquire relevant data of the incidence of COVID-19 at regional level, we have not found reliable data organized by provinces or similar territories so only analysis of the whole country has been made. France is special as it is the reference on which we will base the value for the factor k from Equations (46) and (47) for all countries. We have cautiously chosen a k which is as close as 1 as possible but still generates values for the transmission rate ${\beta }_E$ coherent with the values expected from the initial lockdown stage. In the $\mu$-SEIR model, after the initial lockdown (Stage 1) with ${\beta }_{I_0}={\beta }_{E_0}=0.065$ and the low cases of summer, in which the transmission rates are not calculated, the ${\beta }_E$ in the new normality (Stage 4) rises up to 2–3 times the initial value ${\beta }_E\in$ (0.15–0.20). Stage 3, corresponding to the level of indoor contact, does not seem to change significantly, although it is higher than the following lockdown (Stage 1) starting late October $({\beta }_E\approx$ 0.015–0.100) and the outdoor lockdown (Stage 2) at Christmas season $({\beta }_E\approx$ 0.045–0.130). The variation of the rates in the extended SEIR model is similar as in the $\mu$-SEIR model, although the values of the transmission rates are higher. Here, the initial ${\beta }_{I_0}={\beta }_{E_0}\approx 0.095$ and in the following Stages 4 and 3 ${\beta }_E\in$ (0.265–0.550). During the second lockdown, the calculated ${\beta }_E$ drops to values close to zero when ${\beta }_I={\beta }_{I_0}$, and rises up to ${\beta }_E=0.38$ when ${\beta }_I=0.6{\beta }_{I_0}$. The transmission rate at the final outdoor lockdown Stage (Level 2) presents a infectivity rate ${\beta }_E$ for the exposed quite similar to the previous September-October indoor lockdown Stage (Level 3). In Portugal the variation between the values of ${\beta }_E$ after and before summer is notable. While in the first half of 2020 the values of transmission rates ${\beta }_{E_0},{\beta }_{I_0}$ at the first total lockdown stage in the $\mu$-SEIR and the extended SEIR model are ${\beta }_E=0.062/0.078$ respectively, the first outdoor lockdown (Stage 2) rises up to ${\beta }_E=0.068/0.102$ and in the indoor lockdown (Stage 3) rises again to ${\beta }_E=0.122/0.331$ respectively. The second wave of infection after summer, when the daily cases grows up to 20 daily cases per 100,000 inhabitants, present a ${\beta }_E$ especially higher than the previous one. The values of the transmission rate in the new normality (Stage 4), which extends from the 15th of August to the 15th of September are 0.183 and 0.461 for the $\mu$-SEIR and extended SEIR models respectively; at the following indoor lockdown (Stage 3) they decrease to 0.164 and 0.456, which are values greater than the previous Stage 3 of summer. The second total lockdown, starting the 15th of January 2021 after the rise of new cases in Christmas presents also a value which is higher than the initial total lockdown of 2020, from ${\beta }_E=0.062/0.078$ to ${\beta }_E=0.11/0.42$ for the $\mu$ and extended SEIR models respectively.

6.2. Norway and Sweden

Norway and Sweden present very similar ethnic and sociodemographic profiles and age distributions of the population, health care, educational and political systems [63], population density, age and urban distribution and similar climate due to proximity [64]. However, the impact of the two different strategies applied in both countries is notable as seen in the calculations of the transmission rates ${\beta }_E$ and ${\beta }_I$. Norway proceeded with cautionary NPI measures, opting for the total lockdown at the beginning of the outbreak during March, April and May. After the summer, the rate of new cases per 100,000 inhabitants rose again in both countries and Sweden finally chose as strong of isolation strategy as Norway. The values for the Stage 1, 4 and 3 are, respectively, ${\beta }_E=(0.075,0.530,0.317)/(0.059,0.211,0.114)$ for the $\mu$/extended SEIR model respectively. Sweden, on the other side, while recommending social distancing decided to approach to the situation by appealing to the personal responsibility of the individuals without any NPI measure implemented at the beginning of the outbreak. It can be seen that the increments of new cases in both the two Scandinavian countries (see Figure 4 and Figure 6) appear at the same dates, but the incidence rate of new infections is almost five times higher in Sweden. The transmission rate values for the exposed subpopulation of the new normality Stage 4 and the indoor lockdown (Stage 3) which was practiced after summer in both countries up to February 2021 are quite the same, from ${\beta }_E=0.211/0.113$ as we can see in Figure 5 to ${\beta }_E=0.530/0.317$, in Figure 4, for the $\mu$ and extended SEIR models for Norway, and ${\beta }_E=0.221/0.122$ as we can see in Figure 7 and ${\beta }_E=0.577/0.384$, in Figure 6, in the $\mu$ and extended SEIR models for Sweden. This similarity on this interval suggests that the number of infected individuals up to September, which was approximately 85,000 for a population of 10 million in Sweden in contrast to the 11,000 cases of Norway for a population of 5 million, could have been substantially less in the case of Sweden had they adopted proper NPI measures.

6.3. Spain and Italy

Spain and Italy are the European countries in which COVID-19 initially advanced in 2020 and the countries which acted on the propagation of the disease through NPI earlier in 2020. The fact that the data of incidence of the disease is distributed geographically on the provinces of both countries allows us to make a study at two different levels, a global one and a regional one. The transmission rates ${\beta }_{E_0}$ of Italy and Spain are comparable at the initial total lockdown, which started in March and ended in May as we can see in Figure 18 and Figure 22. In the initial Stage 1, both the infectious and exposed transmission rates are ${\beta }_{I_0}={\beta }_{E_0}=0.065/0.083$ in Italy and ${\beta }_{I_0}={\beta }_{E_0}=0.061/0.078$ in Spain for the $\mu$/extended SEIR model respectively. Observe that in the new normality, after the summer also shows a similar ${\beta }_E=0.268/0.584$ in Italy and ${\beta }_E=0.215/0.543$ in Spain for the $\mu$/extended SEIR model respectively. After the rise of new cases in both countries, a new indoor lockdown (Stage 3) is established in Spain the 15th of August and the 27th of August in Italy. The transmission rates during these stages are reduced to ${\beta }_E=0.148/0.411$ in Spain and ${\beta }_E=0.117/0.347$ in Italy for the $\mu$/extended SEIR models. We can see that while the Stage 3 lockdown NPI were a relative success in Italy as the rate of new cases decreased more rapidly while it kept growing in Spain, this made Spain react earlier in November by reducing down to a Stage 2 lockdown, decreasing the value of the transmission rate to ${\beta }_E=0.04/0.08$ in the $\mu$/extended SEIR model, and the incidence of new infected individuals eventually went down. Italy, on the other hand, maintained the lockdown Stage 3 until the Christmas season, when the decision to go to total lockdown was made, although the period of time in which this measure is applied was not long enough to ${\beta }_E$ to decrease to the levels of the previous one. It can be argued that the cold weather and the special social circumstances of the period in which it was applied affected the probability of infection and thus the ${\beta }_E=0.103/0.167$ at the $\mu$/extended SEIR model. The study of the transmission rates on the 107 and 50 provinces of Italy and Spain show the impact of the different NPI strategies on both countries. On one hand the coherence of the values for ${\beta }_E$ of the first Stage 1 total lockdown in both models and both countries suggests that the measures were more strict at the beginning of the outbreak rather than in the following waves of contagion. On the other hand, the p-test applied to the histograms of the consecutive transition of stages, as seen in the $\mu$/extended SEIR model in Italy (see Table 15 and Table 16) and Spain (see Table 12 and Table 13), confirms that the difference between the values of ${\beta }_E$ is statistically significant and that the impact of the NPIs is real and measurable. The variance of the transmission rates is higher in the new normality and decreases as does the stages of the lockdown.

6.4. Germany and the UK

Germany and the United Kingdom have their data distributed geographically, so the analysis can be made in their regional and country level. While in Germany the implementation of the NPI measures when the second wave of contagion arrived were gradual; in the UK the lockdown measures level changed directly from Level 3 to Level 1. Germany has a transmission rate ${\beta }_E$ 0.060 /0.079, and as the incidence of new cases increases at the beginning of the second wave it rapidly changes to a lockdown Stage 3, where the ${\beta }_E=0.218/0.560$ in the $\mu$/extended SEIR model respectively. However, the rate of new cases is not curved down until the outdoor lockdown (Stage 2) is implemented, when ${\beta }_E=0.134/0.215$ and in the second lockdown at Christmas season ${\beta }_E=0.1/0.215$ in the $\mu$/extended SEIR model. In the following lockdown after Christmas the transmission rate decreases in both the $\mu$/extended model to ${\beta }_E=0.054/0.067$. On the study of the 33 different regions of Germany it can be seen that the variance of the average values of the transmission rates for each is higher during the second lockdown, in contrast with the small variance that occurs at the first total lockdown. As it is seen in Table 9 and Table 10, although the effect of the different isolation measures affects the values of ${\beta }_E$ in a statistical significance, as the p-test made on the histograms shows, the variance of the calculated ${\beta }_E$ during the second wave of the disease increases as the isolation applied is higher. This development, in which the variance in the values of the ${\beta }_E$ of provinces is smaller during the first total lockdown, happens too in the United Kingdom, as it can be seen in the transmission rates of their 64 counties in Table 6 and Table 7 and their related histograms, in which the initial total lockdown has a lower mean value and variance than the following two total lockdowns. The graphics of ${\beta }_E$ of all the UK (see Figure 10 and Figure 12) actually show the impact of these total lockdown stages on the incidence of newly infected individuals, as their numbers decrease, to half in the first transition from indoor lockdown Stage 3 (from August to November) with ${\beta }_E=0.161/0.447$ to total lockdown (November-December) ${\beta }_E=0.078/0.207$ for the $\mu$/extended-SEIR model respectively; and to 1/3 approximately in the second transition from the indoor lockdown Stage 3 (December-January) with ${\beta }_E=0.184/0.485$ to the second total lockdown (since January) with ${\beta }_E=0.49/0.150$.

6.5. Disparities on the Histograms

The first notable issue that comes up when analyzing the results of Spain, Italy, Germany and the United Kingdom and the respective histograms is the differences observed in some phases between the average value of the transmission rates of the local regions and the whole countries. This difference is due to the decrease of the importance of the data from provinces that have greater population. Those provinces add a superior number of cases to the data of the whole country than other less populated areas and, when considering the total country, their importance is higher. In contrast, when making a histogram, the importance of all the provinces is the same. A good example of this problem can be seen in the United Kingdom as more than 17% of the cases from the United Kingdom take place in London and, consequently, a big percentage of all cases are registered in London. The importance of the data from London on the results of the whole United Kingdom is diminished when making an average from the data calculated from the counties, because all of them have the same weight. A weighted arithmetic mean is calculated for the regional rates for each lockdown stage and each country, being the weight for each province proportional to the cases of infection that took place during that stage. Thus, the effect of each province on the average value of the transmission rates ${\beta }_I$ and ${\beta }_E$ depends on the number of cases that have occurred during the phase. We will call the average value the one calculated by the unweighted mean. In Italy, the biggest difference is observed in the value from the total lockdown Stage, i.e., the value of ${\beta }_{I_0}={\beta }_{E_0}$. In Germany, there are not significant differences between the average, the weighted arithmetic mean and the national value of ${\beta }_E$ in all the stages, as all of the states presented on the histogram are big enough to dismiss the disparities. In Spain, as in Italy, the weighted arithmetic mean of ${\beta }_E$ is closer to the value of the whole country in almost all the phases than the average value. The big variance of the average values in Spain, on the other hand, can also be explained as their politics on NPI and isolation measures have been decided on regional levels from early in 2020, and the transition between stages of lockdown has not been coherent through the entire second wave of COVID-19.

7. Conclusions

After analysis of the results, the general perception is that the results are coherent in the sense that the decrements and increments of the infectious rate of the disease are common through all the stages in both models’ studies and in every country. The initial total lockdown of the first wave presents a set of similar values of ${\beta }_{I_0}={\beta }_{E_0}$ around 0.062/0.079 for the $\mu$ and extended SEIR model respectively, with a consistency not only between countries but within the studied provinces. The authors believe this fact is associated to the strong NPI measures that were taken at the beginning of the outbreak due to the lack of knowledge in regards to the disease and the initial alarm this created among society and authorities. The variance of the values of ${\beta }_E$ during the second wave was notably higher than in the first one because of the absence of the same motives as before. The NPI measures got more flexible the moment the authorities were able to cope with the disease without serious risk of sanitary collapse, and thus, the value of ${\beta }_E$ changed through phases more irregularly. Even a study on the same region that applies the same isolation measures would show disparities between the values of the transmission rates in a lockdown stage calculated at different intervals of time, as the methods for measuring the COVID-19 cases, the weather and the attitude about the disease have changed over the seasons of 2020. This suggests that:

1.

A total lockdown really works as a palliative measure in case of an outbreak, as it can be seen when comparing the cases of Sweden and Norway.

2.

The difference in the quality of gathered data, due to availability of tests, management of the health care system and public awareness is present not only between the regions of Europe but also between moments in time.

3.

When no immunization is available, such that big enough vaccination campaigns can be made, any type of NPI measures is able to reduce the contagion in a perceptible way.

The general impression we get from our calculations is that our models work fine enough to describe the situation of the population in regards to the COVID-19 spread. The extended SEIR model on one hand describes all subpopulations in a more accurate way taking into account more precisely the medical observations of the evolution of the disease on an individual, depicting realistically the exposed subpopulation, while the $\mu$-SEIR model is not so accurate, but the model is simpler to operate. It is worth pointing out that although the output rates differ significantly, their relevance is in their relative value, and their changes within lockdown stages vary coherently between each other. The authors hope that future studies on COVID-19 will improve their predictability as a good infrastructure for testing and diagnosis is established, and a better classification of the isolation measures is made. The reduction in the number of diagnoses of other notable airborne transmission diseases, like influenza [65], suggest that seasonality of the disease is a factor to take into account in the impact of any possible outbreaks in the coming months. Although it could be the strictness of the NPI after and before summer, which would change the number of average contacts per individual per day, and thus the transmission ratio of exposed subpopulation, it may also be the influence of the weather on the probability of a contagious interaction between a susceptible and a infected individual [66,67]. The authors believe it to be a mixture of both reasons, which should be studied accordingly. Both models here presented are easily adaptable too when the vaccination campaigns scheduled for 2021 start to have a beneficial impact on the incidence of new cases.