Next Article in Journal
Spatial Fuzzy C-Means Clustering Analysis of U.S. Presidential Election and COVID-19 Related Factors in the Rustbelt States in 2020
Next Article in Special Issue
Design of Type-3 Fuzzy Systems and Ensemble Neural Networks for COVID-19 Time Series Prediction Using a Firefly Algorithm
Previous Article in Journal
Arctan-Based Family of Distributions: Properties, Survival Regression, Bayesian Analysis and Applications
Previous Article in Special Issue
Fractional Modelling and Optimal Control of COVID-19 Transmission in Portugal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Stability Analysis of Delayed COVID-19 Models

by
Mohamed A. Zaitri
,
Cristiana J. Silva
*,† and
Delfim F. M. Torres
Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Axioms 2022, 11(8), 400; https://doi.org/10.3390/axioms11080400
Submission received: 15 July 2022 / Revised: 9 August 2022 / Accepted: 11 August 2022 / Published: 13 August 2022
(This article belongs to the Special Issue Mathematics of the COVID-19)

Abstract

:
We analyze mathematical models for COVID-19 with discrete time delays and vaccination. Sufficient conditions for the local stability of the endemic and disease-free equilibrium points are proved for any positive time delay. The stability results are illustrated through numerical simulations performed in MATLAB.

1. Introduction

A pandemic is an epidemic of an infectious disease occurring worldwide, or over a very wide area, crossing international boundaries and usually affecting a large number of people [1]. An infectious disease, also known as transmissible disease or communicable disease, results from an infection caused by a large range of pathogens such as, for example, bacteria and viruses. Several modes of transmission can be identified, such as droplet, fecal, sexual and/or oral and vector-borne transmissions. Historically, communicable diseases have killed millions of people around the world, for example, smallpox, plague, great flu, polio, tuberculosis, HIV/AIDS, SARS, global H1N1 flu, cholera, measles, and ebola. Following the World Health Organization (WHO), some of these infectious diseases remain a threat for public health, such as cholera, measles, HIV, and tuberculosis [2,3,4].
In December 2019, a very dangerous SARS-CoV-2 virus quickly invaded the city of Wuhan in China and, subsequently, 183 countries in the world [5,6]. WHO declared, on 30 January 2020, the COVID-19 infectious disease as a pandemic [7]. It is now known that the spread of COVID-19 changes very rapidly; therefore, taking appropriate and timely actions can influence the course of this pandemic.
Since the beginning of the COVID-19 pandemic, researchers proposed different and complementary mathematical models that describe, approximately, the spread of SARS-CoV-2 in different regions of the world and with alternative modeling techniques; see, e.g., [8,9,10,11,12,13,14]. Although the literature dealing with models of COVID-19 is now huge, with special issues [15], books [16], and review papers [17] on this topic, deterministic models of COVID-19 with delay differential equations and vaccination are relatively scarce [18].
The introduction of time delays to mathematical epidemic models has been studied in order to better understand and describe the transmission dynamics of infectious diseases; see, e.g., [14,19,20,21]. Moreover, time delays may have an important effect on the stability of the equilibrium points, leading, for example, to periodic solutions by Hopf bifurcation; see, e.g., [22] and references cited therein.
As in other infectious diseases, the latent and incubation periods have an important role in the spread of COVID-19. The latent period of an infectious disease is the time interval between infection and becoming infectious, whether the incubation period is the time interval between infection and the appearance of clinical symptoms [23,24,25]. Following the WHO, the incubation period for COVID-19 is between 2 and 10 days [26]. In [25], the authors estimated the mean latent period to be 5.5 (95% CI: 5.1–5.9) days, shorter than the mean incubation period (6.9 days). However, and differently from other infectious diseases, asymptomatic infected individuals can transmit the infection and this imposes more strict mitigation strategies; see, e.g., [27]. To describe and analyze this biological phenomenon, we generalize here a compartmental mathematical model, first proposed in [28], by considering a system of delayed differential equations with discrete time delays.
In recent years, several epidemic models have been presented, both stochastic and deterministic ones; see e.g., [29,30]. In [28], a deterministic mathematical model is proposed to analyze the spread of the COVID-19 epidemic, based on a dynamic mechanism that incorporates the intrinsic impact of hidden latent and infectious cases on the entire process of virus transmission. In [31], Zaitri et al. applied optimal control theory to a generalized SEIR-type model, based on [28], with three controls, representing social distancing, preventive means, and treatment measures to combat the spread of the COVID-19 pandemic. They analyzed such optimal control problem with respect to real data transmission in Italy. Their results show the appropriateness of the model, in particular with respect to the number of quarantined/hospitalized (confirmed and infected) and recovered individuals. Alternative approaches based on SIR-type models but that combine machine learning methods have also been developed; see, e.g., [32,33].
In our paper, we modify the model analyzed in [28] in order to consider time delays, birth and death rates. More precisely, we introduce a time delay that represents, mathematically, the fact that the migration of individuals from susceptible to infected is subject to delay. Secondly, we present a normalized version of the SEIR-type model, compute the equilibrium points, the basic reproduction number, and we prove sufficient conditions for the stability of the equilibrium points, for any positive time delay. Then, we extend the previous model in order to consider vaccination and perform numerical simulations taking into account the real data of the spread of COVID-19 in Italy from 18 October 2020, to 17 January 2021. This allows us to compare our results with previous ones.
The paper is organized as follows: In Section 2, we propose a delayed S E I Q R P mathematical model for COVID-19. Considering the normalized model of the delayed S E I Q R P model, we prove sufficient conditions for the stability of the equilibrium points for any time delay. Then, in Section 3, we propose a delayed mathematical model for COVID-19 with vaccination. Analogously, we prove sufficient conditions for the stability of the equilibrium points of the normalized s e i q r p w with vaccination, for any time delay. Numerical simulations and a discussion of the results are provided in Section 4, illustrating the stability of both delayed models and their practical utility.

2. The Delayed SEIQRP Model

In this section, we propose a delayed mathematical model for COVID-19, which generalizes the one proposed in [28]. As mentioned in Section 1, there are many different models but, all of them, are approximations of the reality. For example, in [34] the possibility to become susceptible again is ignored, although we know re-infection is possible and occurs; while in [35] deaths are not taken into account.
Our model considers six state variables: susceptible individuals, S ( t ) ; exposed individuals, E ( t ) ; infected individuals, I ( t ) ; quarantined individuals, Q ( t ) ; recovered individuals, R ( t ) ; and insusceptible/protected individuals, P ( t ) . The total population is denoted by N ( t ) and is given by
N ( t ) = S ( t ) + E ( t ) + I ( t ) + Q ( t ) + R ( t ) + P ( t ) , for all t [ 0 , T ] .
The following assumptions are made to describe the spread of COVID-19: b is the birth rate, μ is the death rate, α is the protection rate, β the infection rate, γ the inverse of the average latent time, δ the rate at which infectious people enter in quarantine, and λ the recovery rate. The time delay τ 0 represents the incubation period, that is, the length of time before the infected individuals become infectious.
We introduce a discrete time delay that represents the transfer delay from the class of susceptible individuals to the class of infected individuals, after the contact of a susceptible individual with an infectious one. Precisely, the model we propose is given by the following system of six nonlinear ordinary delayed differential equations:
S ˙ ( t ) = b N ( t ) β S ( t τ ) I ( t τ ) N ( t ) ( α + μ ) S ( t ) , E ˙ ( t ) = β S ( t τ ) I ( t τ ) N ( t ) ( γ + μ ) E ( t ) , I ˙ ( t ) = γ E ( t ) ( δ + μ ) I ( t ) , Q ˙ ( t ) = δ I ( t ) ( λ + μ ) Q ( t ) , R ˙ ( t ) = λ Q ( t ) μ R ( t ) , P ˙ ( t ) = α S ( t ) μ P ( t ) ,
where the state variables are subject to the initial conditions S ( θ ) = S 0 , θ [ τ , 0 ] , E ( 0 ) = E 0 , I ( θ ) = I 0 , θ [ τ , 0 ] , Q ( 0 ) = Q 0 , R ( 0 ) = R 0 , and P ( 0 ) = P 0 .

2.1. The Normalized s e i q r p Delayed Model

In the situation where the total population size N ( t ) is not constant along time, it is often convenient to consider the proportions of each compartment of individuals in the population, namely s ( t ) = S ( t ) N ( t ) , e ( t ) = E ( t ) N ( t ) , i ( t ) = I ( t ) N ( t ) , q ( t ) = Q ( t ) N ( t ) , r ( t ) = R ( t ) N ( t ) , and p ( t ) = P ( t ) N ( t ) . According to equality (1), we have N ˙ ( t ) = ( b μ ) N ( t ) . Therefore, the normalized s e i q r p delayed model is given by
s ˙ ( t ) = b β s ( t τ ) i ( t τ ) ( α + b ) s ( t ) , e ˙ ( t ) = β s ( t τ ) i ( t τ ) ( γ + b ) e ( t ) , i ˙ ( t ) = γ e ( t ) ( δ + b ) i ( t ) , q ˙ ( t ) = δ i ( t ) ( λ + b ) q ( t ) , r ˙ ( t ) = λ q ( t ) b r ( t ) , p ˙ ( t ) = α s ( t ) b p ( t ) .
The state variables for system (3) are subject to the following initial conditions: s ( θ ) = S 0 N ( 0 ) , θ [ τ , 0 ] , e ( 0 ) = E 0 N ( 0 ) , i ( θ ) = I 0 N ( 0 ) , θ [ τ , 0 ] , q ( 0 ) = Q 0 N ( 0 ) , r ( 0 ) = R 0 N ( 0 ) , and p ( 0 ) = P 0 N ( 0 ) , with s ( t ) + e ( t ) + i ( t ) + q ( t ) + r ( t ) + p ( t ) = 1 .
In Section 2.2 we show that model (3) has two equilibrium points: the disease free and the endemic equilibrium.

2.2. Equilibrium Points and the Basic Reproduction Number

The disease free equilibrium and the endemic equilibrium point are obtained by solving the right hand side of equations in (3) equal to zero:
b β s ( t τ ) i ( t τ ) ( α + b ) s ( t ) = 0 , β s ( t τ ) i ( t τ ) ( γ + b ) e ( t ) = 0 , γ e ( t ) ( δ + b ) i ( t ) = 0 , δ i ( t ) ( λ + b ) q ( t ) = 0 , λ q ( t ) b r ( t ) = 0 , α s ( t ) b p ( t ) = 0 ,
from which the disease free equilibrium, Σ 0 , is given by
Σ 0 = s 0 , e 0 , i 0 , q 0 , r 0 , p 0 = b α + b , 0 , 0 , 0 , 0 , α α + b ,
while the endemic equilibrium point, Σ + , is given by
Σ + = s + , e + , i + , q + , r + , p +
with
s + = δ + b γ + b β γ , e + = β s + i + γ + b , i + = β γ b δ + b γ + b α + b β δ + b γ + b , q + = β γ b δ δ δ + b γ + b α + b β λ + b δ + b γ + b , r + = λ δ β γ b λ δ δ + b γ + b α + b b β λ + b δ + b γ + b , p + = α δ + b γ + b b β γ .
Following the method of van den Driessche [36], one easily compute the following basic reproduction number:
R 0 = β γ b α + b δ + b γ + b .
The reader interested in the details of the algorithm according to which the basic reproduction number (7) is computed, is referred to the open access article [37].

2.3. Stability of the Normalized s e i q r p Delayed Model

Now, we prove some sufficient conditions for the local asymptotic stability of the disease free equilibrium, Σ 0 , and the endemic equilibrium point, Σ + , for any time delay τ 0 .
Consider the following coordinate transformation: x 1 ( t ) = s ( t ) s ¯ , x 2 ( t ) = e ( t ) e ¯ , x 3 ( t ) = i ( t ) i ¯ , x 4 ( t ) = q ( t ) q ¯ , x 5 ( t ) = r ( t ) r ¯ , and x 6 ( t ) = p ( t ) p ¯ , where ( s ¯ , r ¯ , i ¯ , q ¯ , r ¯ , p ¯ ) denotes any equilibrium point of system (3). The linearized system of (3) takes the form
X ˙ ( t ) = A 0 X ( t ) + A 1 X ( t τ ) ,
where X = ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) T ,
A 0 = α b 0 0 0 0 0 0 γ b 0 0 0 0 0 γ δ b 0 0 0 0 0 δ λ b 0 0 0 0 0 λ b 0 α 0 0 0 0 b ,
and
A 1 = β i ¯ 0 β s ¯ 0 0 0 β i ¯ 0 β s ¯ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .
The characteristic equation of system (3), for any equilibrium point, is given by
Δ ( y ) = | y I d 6 × 6 A 0 A 1 e τ y | .
We are now in a position to prove our first two results.
Theorem 1
(Stability of the disease free equilibrium of system (3)). If R0 < 1, then the disease free equilibrium Σ0 is locally asymptotically stable for any time-delay τ ≥ 0. If R0 > 1, then the disease free equilibrium Σ0 is unstable for any time-delay τ ≥ 0.
Proof. 
The characteristic equation of (3), at the disease free equilibrium Σ 0 , is given by
P ( y , τ ) = ( y + b ) 2 ( y + α + b ) ( y + λ + b ) ( y 2 + Λ 1 y + Λ 2 ( y ) ) = 0 ,
where Λ 1 = δ + 2 b + γ and Λ 2 ( y ) = ( δ + b ) ( γ + b ) β γ b α + b e τ y .
Let R 0 < 1 . We divide the proof into the non-delayed and delayed cases.
(i)
Let τ = 0 . In this case, the Equation (10) becomes
P ( y , 0 ) = ( y + b ) 2 ( y + α + b ) ( y + λ + b ) y 2 + Λ 1 y + ( δ + b ) ( γ + b ) β γ b α + b = 0 .
We need to prove that all roots of the characteristic Equation (11) have negative real parts. It is easy to see that y 1 = b , y 2 = α b and y 3 = λ b are roots of Equation (11) and all of them are real negative roots. Thus, we just need to analyze the fourth term of (11), here denoted by P 1 , that is,
P 1 ( y , 0 ) : = y 2 + Λ 1 y + ( δ + b ) ( γ + b ) β γ b α + b .
Using the Routh–Hurwitz criterion [38], we know that all roots of P 1 ( y , 0 ) have negative real parts if, and only if, the coefficients of P 1 ( y , 0 ) are strictly positive. In this case, we have Λ 1 = δ + 2 b + γ > 0 and
( δ + b ) ( γ + b ) β γ b α + b > 0 if   and   only   if R 0 = β γ b δ + b γ + b α + b < 1 .
Therefore, we have just proved that the disease free equilibrium, Σ 0 , is locally asymptotically stable for τ = 0 , whenever R 0 < 1 .
(ii)
Let τ > 0 . In this case, we will use Rouché’s theorem [39,40] to prove that all roots of the characteristic Equation (10) cannot intersect the imaginary axis, i.e., the characteristic equation cannot have pure imaginary roots. Suppose the contrary, that is, suppose there exists w R such that y = w i is a solution of (10). Replacing y in the fourth term of (10), we get that
w 2 + ( δ + 2 b + γ ) w i + ( δ + b ) ( γ + b ) β γ b α + b cos ( τ w ) i sin ( τ w ) = 0 .
Then,
w 2 + ( δ + b ) ( γ + b ) = β γ b α + b cos ( τ w ) , ( δ + 2 b + γ ) w = β γ b α + b sin ( τ w ) .
By adding up the squares of both equations, and using the fundamental trigonometric formula, we obtain that
w 4 + δ + b 2 + γ + b 2 w 2 + ( δ + b ) 2 ( γ + b ) 2 β γ b α + b 2 = 0 ,
which is equivalent to
w 2 = 1 2 δ + b 2 γ + b 2 2 + 4 β γ b α + b 2 1 2 δ + b 2 + γ + b 2 .
If R 0 < 1 , then ( δ + b ) 2 ( γ + b ) 2 β γ b α + b 2 > 0 , and
δ + b 2 + γ + b 2 2 4 ( δ + b ) 2 ( γ + b ) 2 β γ b α + b 2 < δ + b 2 + γ + b 2 2 ,
so that
δ + b 2 γ + b 2 2 + 4 β γ b α + b 2 < δ + b 2 + γ + b 2 .
Hence, we have w 2 < 0 , which is a contradiction. Therefore, we have proved that whenever R 0 < 1 , the characteristic Equation (10) cannot have pure imaginary roots and the disease free equilibrium Σ 0 is locally asymptotically stable, for any strictly positive time-delay τ .
(iii)
Suppose now that R 0 > 1 . We know that the characteristic Equation (10) has three real negative roots y 1 = b , y 2 = α b , and y 3 = λ b . Thus, we need to check if the remaining roots of
q ( y ) : = y 2 + Λ 1 y + Λ 2 ( y )
have negative real parts. It is easy to see that q ( 0 ) = Λ 2 ( 0 ) < 0 because we are assuming R 0 > 1 . On the other hand, lim y + q ( y ) = + . Therefore, by continuity of q ( y ) , there is at least one positive root of the characteristic Equation (10). Hence, we conclude that Σ 0 is unstable when R 0 > 1 .
The proof is complete. □
Theorem 2
(Stability of the endemic equilibrium point of system (3)). Let τ = 0 . If R 0 > 1 , then the endemic equilibrium point Σ + is locally asymptotically stable. When τ > 0 , the endemic equilibrium point Σ + is locally asymptotically stable if the basic reproduction number R 0 satisfies the following relations:
1 < R 0 < min 3 , 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b
and
M 1 R 0 2 + M 2 R 0 + M 3 > 0 ,
where
M 1 = ( α + b ) 2 ( δ + b ) 2 + ( γ + b ) 2 , M 2 = 2 α + b 2 δ + γ + 2 b 2 3 ( δ + b ) ( γ + b ) + 2 α + b ( δ + b ) ( γ + b ) ( δ + γ + 2 b ) , M 3 = 2 ( α + b ) ( δ + b ) ( γ + b ) α δ γ b .
Proof. 
The characteristic Equation (9), computed at the endemic equilibrium point Σ + , is given by
P ˜ ( y , τ ) = ( y + b ) 2 ( y + λ + b ) ( y 3 + Δ 1 ( y ) y 2 + Δ 2 ( y ) y + Δ 3 ( y ) ) = 0 ,
where Δ 1 ( y ) = L 1 + L ¯ 1 e τ y , Δ 2 ( y ) = L 2 + L ¯ 2 e τ y , and Δ 3 ( y ) = L 3 + L ¯ 3 e τ y with
L 1 = α + δ + γ + 3 b , L ¯ 1 = β γ b δ + b γ + b α + b δ + b γ + b , L 2 = ( γ + 2 b + δ ) ( α + b ) + ( γ + b ) ( δ + b ) , L ¯ 2 = ( γ + 2 b + δ ) ( α + b ) ( R 0 1 ) ( γ + b ) ( δ + b ) , L 3 = ( α + b ) ( γ + b ) ( δ + b ) , L ¯ 3 = β γ b 2 ( α + b ) ( γ + b ) ( δ + b ) .
(i)
Let τ = 0 . In this case, the Equation (16) becomes
P ˜ ( y , 0 ) = ( y + b ) 2 ( y + λ + b ) ( y 3 + Δ ˜ 1 y 2 + Δ ˜ 2 y + Δ ˜ 3 ) = 0 ,
where Δ ˜ 1 = L 1 + L ¯ 1 , Δ ˜ 2 = L 2 + L ¯ 2 and Δ ˜ 3 = L 3 + L ¯ 3 . We need to prove that all the roots of the characteristic Equation (17) have negative real parts. It is easy to see that y 1 = b and y 2 = λ b are roots of (17) and both are real negative roots. Thus, we just need to consider the third term of the above equation. Let
P ˜ 3 ( y , 0 ) : = y 3 + Δ ˜ 1 y 2 + Δ ˜ 2 y + Δ ˜ 3 = 0 .
Using the Routh–Hurwitz criterion [38], we know that all roots of P ˜ 3 ( y , 0 ) have negative real parts if, and only if, the coefficients of P ˜ 3 ( y , 0 ) are strictly positive and Δ ˜ * = Δ ˜ 1 Δ ˜ 2 Δ ˜ 3 > 0 . If R 0 > 1 , then
Δ ˜ 1 = α + δ + γ + 3 b + ( α + b ) ( R 0 1 ) > 0 , Δ ˜ 2 = ( δ + γ + 2 b ) ( α + b ) R 0 > 0 , Δ ˜ 3 = ( α + b ) ( δ + b ) ( γ + b ) ( R 0 1 ) > 0 , Δ ˜ * = ( α + b ) ( α + b ) ( δ + γ + 2 b ) R 0 2 + ( α + b ) ( δ 2 + 3 b ( δ + b ) + γ ( δ + γ + 3 b ) ) R 0 + ( α + b ) ( δ + b ) ( γ + b ) > 0 .
(ii)
Let τ > 0 . Using Rouché’s theorem, we prove that all the roots of the characteristic Equation (16) cannot intersect the imaginary axis, i.e., the characteristic equation cannot have pure imaginary roots. Suppose the opposite, that is, assume there exists w R such that y = w i is a solution of (16). Replacing y into the third term of (16), we get that
w 3 i L 1 w 2 + L 2 w i + L 3 + ( L ¯ 1 w 2 + L ¯ 2 w i + L ¯ 3 ) cos ( τ w ) i sin ( τ w ) = 0 .
Then,
L 1 w 2 + L 3 = ( L ¯ 1 w 2 L ¯ 3 ) cos ( τ w ) L ¯ 2 w sin ( τ w ) , w 3 + L 2 w = L ¯ 2 w cos ( τ w ) ( L ¯ 1 w 2 L ¯ 3 ) sin ( τ w ) .
By adding up the squares of both equations, and using the fundamental trigonometric formula, we obtain that
w 6 + K 1 w 4 + K 2 w 2 + K 3 = 0 ,
where
K 1 = L 1 2 L ¯ 1 2 2 L 2 , K 2 = 2 L ¯ 1 L ¯ 3 2 L 1 L 3 + L 2 2 L ¯ 2 2 , K 3 = L 3 2 L ¯ 3 2 .
Assume that the basic reproduction number R 0 satisfies relations (14) and (15) with the following condition:
min 3 , 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b = 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b .
Then,
K 1 = ( δ + b ) 2 + ( γ + b ) 2 + ( α + b ) 2 1 R 0 1 2 > 0 .
In contrast, if R 0 satisfies relations (14) and (15) with the condition
min 3 , 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b = 3 ,
then we have
1 < R 0 < 3 < 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b ,
which is equivalent to
0 < R 0 1 < 2 < ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b , 1 ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 ( α + b ) 2 < 1 ( R 0 1 ) 2 < 1 , ( δ + b ) 2 ( γ + b ) 2 < ( α + b ) 2 1 ( R 0 1 ) 2 < ( α + b ) 2 .
Thus,
K 1 > 0 .
Under the assumption that the basic reproduction number R 0 satisfies relations (14) and (15), we have
K 2 = M 1 R 0 2 + M 2 R 0 + M 3 > 0 .
Therefore, if we assume that the basic reproduction number R 0 satisfies relations (14) and (15) with condition (20), then
K 3 = ( α + b ) 2 ( δ + b ) 2 ( γ + b ) 2 1 R 0 2 2 > 0 ;
if R 0 satisfies relations (14) and (15) with condition (19), then we have
1 < R 0 < 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b < 3 ,
which is equivalent to
1 < R 0 2 < 1 + ( α + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b < 1 ,
and also equivalent to
1 R 0 2 2 > 0 .
Thus,
K 3 > 0 .
We conclude that the left hand-side of equation (16) is strictly positive, which implies that this equation is not possible. Therefore, (17) does not have imaginary roots, which implies that Σ + is locally asymptotically stable for any time delay τ > 0 .
The proof is complete. □
It should be noted that Theorem 2 is not trivial, and it is not easy to give a biological/medical interpretation to the relations (14) and (15).

3. The Delayed SEIQRPW Model with Vaccination

Let us introduce in model (2) a constant u and an extra variable W ( t ) , t [ 0 , t f ] , representing the fraction of susceptible individuals that are vaccinated and the number of vaccines used, respectively, with
W ˙ ( t ) = u S ( t ) ,
subject to the initial condition W ( 0 ) = 0 . Note that (21) is just the production rate of vaccinated.
The model with vaccination is given by the following system of seven nonlinear delayed differential equations:
S ˙ ( t ) = b N ( t ) β S ( t τ ) I ( t τ ) N ( t ) α + u + μ S ( t ) , E ˙ ( t ) = β S ( t τ ) I ( t τ ) N ( t ) ( γ + μ ) E ( t ) , I ˙ ( t ) = γ E ( t ) ( δ + μ ) I ( t ) , Q ˙ ( t ) = δ I ( t ) ( λ + μ ) Q ( t ) , R ˙ ( t ) = λ Q ( t ) μ R ( t ) , P ˙ ( t ) = α S ( t ) μ P ( t ) , W ˙ ( t ) = u S ( t ) μ W ( t ) ,
where the total population N ( t ) is given by
N ( t ) = S ( t ) + E ( t ) + I ( t ) + Q ( t ) + R ( t ) + P ( t ) + W ( t ) , t [ 0 , T ] .
The state variables are subject to the following initial conditions: S ( θ ) = S 0 , θ [ τ , 0 ] , E ( 0 ) = E 0 , I ( θ ) = I 0 , θ [ τ , 0 ] , Q ( 0 ) = Q 0 , R ( 0 ) = R 0 , P ( 0 ) = P 0 , and W ( 0 ) = 0 .
Note that in model (22) we do not vaccinate the insusceptible/protected individuals P ( t ) , assumed protected through precautionary measures with a protection rate α . Moreover, the fraction of susceptible individuals that are vaccinated is u.

3.1. Normalized s e i q r p w Delayed Model with Vaccination

Analogously to Section 2, we consider the proportions of each compartment of individuals in the population, namely s ( t ) = S ( t ) N ( t ) , e ( t ) = E ( t ) N ( t ) , i ( t ) = I ( t ) N ( t ) , q ( t ) = Q ( t ) N ( t ) , r ( t ) = R ( t ) N ( t ) , p ( t ) = P ( t ) N ( t ) , and w ( t ) = W ( t ) N ( t ) . According to Equation (23), we have N ˙ ( t ) = ( b μ ) N ( t ) . Therefore, the normalized s e i q r p w delayed model is given by
s ˙ ( t ) = b β s ( t τ ) i ( t τ ) ( α + u + b ) s ( t ) , e ˙ ( t ) = β s ( t τ ) i ( t τ ) ( γ + b ) e ( t ) , i ˙ ( t ) = γ e ( t ) ( δ + b ) i ( t ) , q ˙ ( t ) = δ i ( t ) ( λ + b ) q ( t ) , r ˙ ( t ) = λ q ( t ) b r ( t ) , p ˙ ( t ) = α s ( t ) b p ( t ) , w ˙ ( t ) = u s ( t ) b w ( t ) .
The state variables for system (24) are subject to the following initial conditions: s ( θ ) = S 0 N ( 0 ) , θ [ τ , 0 ] , e ( 0 ) = E 0 N ( 0 ) , i ( θ ) = I 0 N ( 0 ) , θ [ τ , 0 ] , q ( 0 ) = Q 0 N ( 0 ) , r ( 0 ) = R 0 N ( 0 ) , p ( 0 ) = P 0 N ( 0 ) , and w ( 0 ) = 0 , with s ( t ) + e ( t ) + i ( t ) + q ( t ) + r ( t ) + p ( t ) + w ( t ) = 1 .

3.2. Equilibrium Points and the Basic Reproduction Number

The disease free and the endemic equilibrium points of model (24) can be obtained by equating the right-hand side of Equation (24) to zero, hence satisfying
b β s ( t τ ) i ( t τ ) ( α + u + b ) s ( t ) = 0 , β s ( t τ ) i ( t τ ) ( γ + b ) e ( t ) = 0 , γ e ( t ) ( δ + b ) i ( t ) = 0 , δ i ( t ) ( λ + b ) q ( t ) = 0 , λ q ( t ) b r ( t ) = 0 , α s ( t ) b p ( t ) = 0 , u s ( t ) b w ( t ) = 0 .
The disease free equilibrium of model (24), Σ 1 , is given by
Σ 1 = s 0 , e 0 , i 0 , q 0 , r 0 , p 0 , w 0 = b α + u + b , 0 , 0 , 0 , 0 , α α + u + b , u α + u + b ,
while the endemic equilibrium point for system (24), Σ V + , is given by
Σ V + = s * , e * , i * , q * , r * , p * , w * ,
where
s * = δ + b γ + b β γ , e * = β s + i + γ + b , i * = β γ b δ + b γ + b α + u + b β δ + b γ + b , q * = β γ b δ δ δ + b γ + b α + u + b β λ + b δ + b γ + b , r * = λ δ β γ b λ δ δ + b γ + b α + u + b b β λ + b δ + b γ + b , p * = α δ + b γ + b b β γ , w * = u δ + b γ + b b β γ .
Following the method from van den Driessche [36], we obtain the following basic reproduction number, denoted by R ˜ 0 :
R ˜ 0 = β γ b α + u + b δ + b γ + b .

3.3. Stability of the Normalized s e i q r p w Delayed Model with Vaccination

Consider the following coordinate transformation: x 1 ( t ) = s ( t ) s ¯ , x 2 ( t ) = e ( t ) e ¯ , x 3 ( t ) = i ( t ) i ¯ , x 4 ( t ) = q ( t ) q ¯ , x 5 ( t ) = r ( t ) r ¯ , x 6 ( t ) = p ( t ) p ¯ , and x 7 ( t ) = w ( t ) w ¯ , where ( s ¯ , r ¯ , i ¯ , q ¯ , r ¯ , p ¯ , w ¯ ) denotes an equilibrium point of system (24). The linearized system of (24) takes the form
X ˙ ( t ) = A ˜ 0 X ( t ) + A ˜ 1 X ( t τ ) ,
where X = ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) T ,
A ˜ 0 = α u b 0 0 0 0 0 0 0 γ b 0 0 0 0 0 0 γ δ b 0 0 0 0 0 0 δ λ b 0 0 0 0 0 0 λ b 0 0 α 0 0 0 0 b 0 u 0 0 0 0 0 b ,
A ˜ 1 = β i ¯ 0 β s ¯ 0 0 0 0 β i ¯ 0 β s ¯ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .
The characteristic equation of system (24) is given by
Γ ˜ ( y ) = | y I d 7 × 7 A ˜ 0 A ˜ 1 e τ y | .
We are also able to prove stability results for the normalized s e i q r p w delayed model with vaccination.
Theorem 3
(Stability of the disease free equilibrium of system (24)). If R ˜ 0 < 1 , then the disease free equilibrium Σ 1 is locally asymptotically stable for any time-delay τ 0 . If R ˜ 0 > 1 , then the disease free equilibrium is unstable for any time-delay τ 0 .
Proof. 
The characteristic Equation (29) at the disease free equilibrium, Σ 1 , is given by
P * ( y , τ ) = ( y + b ) 3 ( y + α + u + b ) ( y + λ + b ) ( y 2 + Γ 1 y + Γ 2 ( y ) ) = 0 ,
where Γ 1 = δ + 2 b + γ and Γ 2 ( y ) = ( δ + b ) ( γ + b ) β γ b α + u + b e τ y .
(i)
Let τ = 0 . In this case, the Equation (30) becomes
P * ( y , 0 ) = ( y + b ) 3 ( y + α + u + b ) ( y + λ + b ) y 2 + Γ 1 y + ( δ + b ) ( γ + b ) β γ b α + u + b = 0 .
We need to prove that all roots of the characteristic Equation (31) have negative real parts. It is easy to see that y 1 = b , y 2 = α u b and y 3 = λ b are roots of Equation (31) and the three are real and negative. Thus, we just need to consider the fourth term of Equation (31). Let
P 3 * ( y , 0 ) : = y 2 + Γ 1 y + ( δ + b ) ( γ + b ) β γ b α + u + b .
Using the Routh–Hurwitz criterion [38], we know that all roots of P 3 * ( y , 0 ) have negative real parts if, and only if, the coefficients of P 3 * ( y , 0 ) are strictly positive. In this case, Γ 1 = δ + 2 b + γ > 0 and
( δ + b ) ( γ + b ) β γ b α + u + b > 0 if , and only if , R ˜ 0 = β γ b α + u + b δ + b γ + b < 1 .
Therefore, we have proved that the disease free equilibrium, Σ 1 , is locally asymptotically stable for τ = 0 , whenever R ˜ 0 < 1 .
(ii)
Let τ > 0 . Using Rouché’s theorem, we prove that all roots of the characteristic Equation (30) cannot have pure imaginary roots. Suppose the contrary, i.e., that there exists w R such that y = w i is a solution of (30). Replacing y in the fourth term of (30), we get
w 2 + ( δ + 2 b + γ ) w i + ( δ + b ) ( γ + b ) β γ b α + u + b cos ( τ w ) i sin ( τ w ) = 0 .
Then,
w 2 + ( δ + b ) ( γ + b ) = β γ b α + u + b cos ( τ w ) , ( δ + 2 b + γ ) w = β γ b α + u + b sin ( τ w ) .
By adding up the squares of both equations and using the fundamental trigonometric formula, one has
w 4 + δ + b 2 + γ + b 2 w 2 + ( δ + b ) 2 ( γ + b ) 2 β γ b α + u + b 2 = 0 ,
which is equivalent to
w 2 = 1 2 δ + b 2 γ + b 2 2 + 4 β γ b α + u + b 2 1 2 δ + b 2 + γ + b 2 .
If R ˜ 0 < 1 , then ( δ + b ) 2 ( γ + b ) 2 β γ b α + u + b 2 > 0 , and
δ + b 2 + γ + b 2 2 4 ( δ + b ) 2 ( γ + b ) 2 β γ b α + u + b 2 < δ + b 2 + γ + b 2 2 ,
so that
δ + b 2 γ + b 2 2 + 4 β γ b α + u + b 2 < δ + b 2 + γ + b 2 .
Hence, we have w 2 < 0 , which is a contradiction. Therefore, we have proved that if R ˜ 0 < 1 , then the characteristic Equation (30) cannot have pure imaginary roots and the disease free equilibrium Σ 1 is locally asymptotically stable, for any strictly positive time delay τ .
(iii)
Suppose now that R ˜ 0 > 1 . We know that the characteristic Equation (30) has three real negative roots y 1 = b , y 2 = α u b and y 3 = λ b . Thus, we need to check if the remaining roots of
q * ( y ) : = y 2 + Γ 1 y + Γ 2 ( y )
have negative real parts. It is easy to see that q ( 0 ) = Γ 2 ( 0 ) < 0 , because we are assuming R ˜ 0 > 1 . On the other hand, lim y + q * ( y ) = + . Therefore, by continuity of q * ( y ) , there is at least one positive root of the characteristic Equation (30). Hence, we conclude that Σ 1 is unstable, for any τ 0 .
The proof is complete. □
Theorem 4
(Stability of the endemic equilibrium point of system (24)). Let τ = 0 . If R ˜ 0 > 1 , then the endemic equilibrium point Σ V + is locally asymptotically stable. When τ > 0 , the endemic equilibrium point Σ V + is locally asymptotically stable if the basic reproduction number R ˜ 0 satisfies the following relations:
1 < R ˜ 0 < min 3 , 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b
and
M 1 * R ˜ 0 2 + M 2 * R ˜ 0 + M 3 * > 0 ,
where
M 1 * = ( α + u + b ) 2 ( δ + b ) 2 + ( γ + b ) 2 , M 2 * = 2 α + u + b 2 δ + γ + 2 b 2 3 ( δ + b ) ( γ + b ) + 2 α + u + b ( δ + b ) ( γ + b ) ( δ + γ + 2 b ) , M 3 * = 2 ( α + u + b ) ( δ + b ) ( γ + b ) α + u δ γ b .
Proof. 
The characteristic Equation (29), computed at the endemic equilibrium Σ V + , is given by
P ˜ * ( y , τ ) = ( y + b ) 3 ( y + λ + b ) ( y 3 + Ω 1 ( y ) y 2 + Ω 2 ( y ) y + Ω 3 ( y ) ) = 0 ,
where Ω 1 ( y ) = L 1 * + L ¯ 1 * e τ y , Ω 2 ( y ) = L 2 * + L ¯ 2 * e τ y , and Ω 3 ( y ) = L 3 * + L ¯ 3 * e τ y with
L 1 * = α + u + δ + γ + 3 b , L ¯ 1 * = β γ b δ + b γ + b α + u + b δ + b γ + b , L 2 * = ( γ + 2 b + δ ) ( α + u + b ) + ( γ + b ) ( δ + b ) , L ¯ 2 * = ( γ + 2 b + δ ) ( α + u + b ) ( R ˜ 0 1 ) ( γ + b ) ( δ + b ) , L 3 * = ( α + u + b ) ( γ + b ) ( δ + b ) , L ¯ 3 * = β γ b 2 ( α + u + b ) ( γ + b ) ( δ + b ) .
(i)
Let τ = 0 . In this case, Equation (37) becomes
P ˜ * ( y , 0 ) = ( y + b ) 3 y + λ + b y 3 + Ω ˜ 1 y 2 + Ω ˜ 2 y + Ω ˜ 3 = 0 ,
where Ω ˜ 1 = L 1 * + L ¯ 1 * , Ω ˜ 2 = L 2 * + L ¯ 2 * and Ω ˜ 3 = L 3 * + L ¯ 3 * . Looking at the roots of the characteristic Equation (38), it is easy to see that y 1 = b and y 2 = λ b are real negative roots of (38). Considering the third term of the above equation, let
P ˜ 3 * ( y , 0 ) : = y 3 + Ω ˜ 1 y 2 + Ω ˜ 2 y + Ω ˜ 3 = 0 .
Using the Routh–Hurwitz criterion [38], we know that all roots of P ˜ 3 * ( y , 0 ) have negative real parts if, and only if, the coefficients of P ˜ 3 * ( y , 0 ) are strictly positive and
Ω ˜ * = Ω ˜ 1 Ω ˜ 2 Ω ˜ 3 > 0 .
If R ˜ 0 > 1 , then
Ω ˜ 1 = α + u + δ + γ + 3 b + α + u + b ( R ˜ 0 1 ) > 0 , Ω ˜ 2 = ( δ + γ + 2 b ) ( α + u + b ) R ˜ 0 > 0 , Ω ˜ 3 = ( α + u + b ) ( δ + b ) ( γ + b ) ( R ˜ 0 1 ) > 0 , Ω ˜ * = ( α + u + b ) ( δ + γ + 2 b ) R ˜ 0 2 + ( α + u + b ) ( δ 2 + 3 b ( δ + b ) + γ ( δ + γ + 3 b ) ) R ˜ 0 + ( α + u + b ) ( δ + b ) ( γ + b ) > 0 .
(ii)
Let τ > 0 . By Rouché’s theorem, we prove that all roots of the characteristic Equation (37) cannot intersect the imaginary axis, i.e., the characteristic equation cannot have pure imaginary roots. Suppose the opposite, i.e., that there exists w R such that y = w i is a solution of (37). Replacing y in the third term of (37), we get
w 3 i L 1 * w 2 + L 2 * w i + L 3 * + ( L ¯ 1 * w 2 + L ¯ 2 * w i + L ¯ 3 * ) cos ( τ w ) i sin ( τ w ) = 0 .
Then,
L 1 * w 2 + L 3 * = ( L ¯ 1 * w 2 L ¯ 3 * ) cos ( τ w ) L ¯ 2 * w sin ( τ w ) , w 3 + L 2 * w = L ¯ 2 * w cos ( τ w ) ( L ¯ 1 * w 2 L ¯ 3 * ) sin ( τ w ) .
By adding up the squares of both equations and using the fundamental trigonometric formula, we obtain that
w 6 + K 1 * w 4 + K 2 * w 2 + K 3 * = 0 ,
where
K 1 * = ( L 1 * ) 2 ( L ¯ 1 * ) 2 2 L 2 * , K 2 * = 2 L ¯ 1 * L ¯ 3 * 2 L 1 * L 3 * + ( L 2 * ) 2 ( L ¯ 2 * ) 2 , K 3 * = ( L 3 * ) 2 ( L ¯ 3 * ) 2 .
Assume that the basic reproduction number R ˜ 0 satisfies relations (34) and (35) with the condition
min 3 , 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + b = 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b .
Then,
K 1 * = ( δ + b ) 2 + ( γ + b ) 2 + ( α + u + b ) 2 ( 1 ( R ˜ 0 1 ) 2 ) > 0 .
In contrast, if R ˜ 0 satisfies relations (34) and (35) under the condition
min 3 , 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b = 3 ,
then we have
1 < R ˜ 0 < 3 < 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b ,
which is equivalent to
0 < R ˜ 0 1 < 2 < ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b , 1 ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 ( α + u + b ) 2 < 1 ( R ˜ 0 1 ) 2 < 1 , ( δ + b ) 2 ( γ + b ) 2 < ( α + u + b ) 2 1 ( R ˜ 0 1 ) 2 < ( α + u + b ) 2 .
Thus,
K 1 * > 0 .
Under the assumption that the basic reproduction number R ˜ 0 satisfies relations (34) and (35), we have
K 2 * = M 1 * R ˜ 0 2 + M 2 * R ˜ 0 + M 3 * > 0 .
Therefore, if we assume that the basic reproduction number R ˜ 0 satisfies relations (34) and (35) with condition (41), then
K 3 * = ( α + u + b ) 2 ( δ + b ) 2 ( γ + b ) 2 1 R ˜ 0 2 2 > 0 ;
if R ˜ 0 satisfies (34) and (35) with condition (40), then we have
1 < R ˜ 0 < 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b < 3 ,
which is equivalent to
1 < R ˜ 0 2 < 1 + ( α + u + b ) 2 + ( δ + b ) 2 + ( γ + b ) 2 α + u + b < 1 ,
and also equivalent to
1 R ˜ 0 2 2 > 0 .
Thus,
K 3 * > 0 .
We have just proved that the left hand-side of Equation (37) is strictly positive, which implies that this equation is not possible. Therefore, (38) does not have imaginary roots, and Σ V + is locally asymptotically stable, for any time delay τ > 0 , whenever R ˜ 0 satisfies conditions (34) and (35).
The proof is complete. □
It should be noted that Theorem 3 is not trivial, and it is not easy to give a biological/medical interpretation to the relations (34) and (35).

4. Numerical Simulations and Discussion

In this section we investigate, numerically, the local stability of the normalized s e i q r p and s e i q r p w models, illustrating our results from Section 2 and Section 3. All numerical computations were performed in the numeric computing environment MATLAB R2019b using the medium order method and numerical interpolation [41].

4.1. Local Stability of the Delayed s e i q r p Model

Consider the normalized delayed s e i q r p model (3), proposed in Section 2. Take the initial conditions
( s 0 , e 0 , i 0 , q 0 , r 0 , p 0 ) = 0.7 , 0.05 , 0.05 , 0.1 , 0.05 , 0.05
and the parameter values as given in Table 1.
In Figure 1, we present the numerical solutions to the delayed model (3) in the time interval [ 0 , 30 ] days.
Considering the parameter values from Table 1, we have the following value for the basic reproduction number R 0 of Section 2: R 0 = 1.5 . From Theorem 2, R 0 = 1.5 satisfies the conditions (14) and (15), so the endemic equilibrium point E E = ( 1 3 , 1 6 , 1 12 , 1 24 , 1 24 , 1 3 ) of system (3) is locally asymptotically stable for any time delay τ 0 .
In Figure 2, we observe the effect of the time delays: τ = 0 , , 6 on the classes e of exposed and i of infectious. The presence of waves is due to the presence of the time delay and is related to the emergence of the COVID-19 waves. For the study of multiple epidemic waves in the context of COVID-19, we refer the interested reader to [42].

4.2. Delayed s e i q r p w Model with Vaccination: COVID-19 in Italy

Now, we study, numerically, the stability of the spread of the epidemic of COVID-19 in Italy for the period of three months starting from 18 October 2020, using the delayed model (24) that we proposed in Section 3. The preliminary conditions and real data were taken and computed from the database https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv (accessed on 14 August 2021). We consider the initial conditions
( s 0 , e 0 , i 0 , q 0 , r 0 , p 0 , w 0 ) = 1 N 59769273 , 403601 , 8837 , 44098 , 254058 , 133 , 0
with N = 60480000 [43], and the parameter values as given in Table 2, which are obtained using the nonlinear least-squares solver [44]. The reader interested in the details of the nonlinear least-squares solver, according to which the parameters of the delayed model (24) are computed, is referred to the open access article [44].
In Figure 3 and Figure 4, we present numerical solutions to the delayed model (24) in the time interval t [ 0 , 90 ] days, t = 0 representing 18 October 2020, and considering two cases:
  • Case 1: τ = 0 days (without delay), with different percentages of susceptible individuals being vaccinated — u = 0 % , u = 20 % , u = 40 % and u = 60 % (Figure 3).
  • Case 2: u = 20 % (fixed), with different delays — τ = 0 days, τ = 3 days, and τ = 6 days (Figure 4).
Considering the parameter values from Table 2, and u = 0 , u = 20 % , u = 40 % , u = 60 % , we have the following values for the basic reproduction number R ˜ 0 of Section 3: R ˜ 0 = 0.0647 , R ˜ 0 = 0.0554 , R ˜ 0 = 0.0484 , and R ˜ 0 = 0.043 , respectively. From Theorem 3, the disease free equilibrium Σ 1 of system (24) is locally asymptotically stable for the time delay τ = 0 . From Theorem 4, the endemic equilibrium point Σ V + of system (24) is unstable for the time delay τ = 0 .
In conclusion, there is an inverse proportional relationship between the fraction u of susceptible individuals that are vaccinated and the number of exposed, infected, and recovered individuals: the greater the fraction of susceptible individuals that are vaccinated, the smaller the number of exposed, infected, and recovered individuals would be, and vice versa (see Figure 3). Moreover, there is a directly proportional relationship between the transfer time delay τ from the class of susceptible individuals to the class of infected individuals and the number of exposed, infected, and recovered individuals: the greater the time delay, the greater the number of exposed, infected, and recovered individuals would be, and vice versa (see Figure 4).

Author Contributions

Conceptualization, M.A.Z., C.J.S. and D.F.M.T.; methodology, M.A.Z., C.J.S. and D.F.M.T.; software, M.A.Z.; validation, M.A.Z., C.J.S. and D.F.M.T.; formal analysis, M.A.Z., C.J.S. and D.F.M.T.; investigation, M.A.Z., C.J.S. and D.F.M.T.; writing—original draft preparation, M.A.Z., C.J.S. and D.F.M.T.; writing—review and editing, M.A.Z., C.J.S. and D.F.M.T.; visualization, M.A.Z.; supervision, C.J.S. and D.F.M.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by FCT (Fundação para a Ciência e a Tecnologia), grant number UIDB/04106/2020 (CIDMA). C.J.S. was also supported by FCT via the FCT Researcher Program CEEC Individual 2018 with reference CEECIND/00564/2018.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. These data can be found here: https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv (accessed on 14 August 2021).

Acknowledgments

The authors are very grateful to four reviewers for several constructive comments, suggestions and questions that helped them to improve their manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Porta, M. A Dictionary of Epidemiology; Oxford University Press: Oxford, UK, 2014. [Google Scholar] [CrossRef]
  2. WHO, World Health Organization. Available online: http://www.emro.who.int/pandemic-epidemic-diseases/outbreaks/index.html (accessed on 29 December 2021).
  3. WHO, World Health Organization. Available online: https://www.who.int/health-topics/hiv-aids (accessed on 29 December 2021).
  4. WHO, World Health Organization. Available online: https://www.who.int/health-topics/tuberculosis (accessed on 29 December 2021).
  5. Lai, C.-C.; Shih, T.P.; Ko, W.C.; Tang, H.J.; Hsueh, P.R. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges. Int. J. Antimicrob. Agents 2020, 55, 105924. [Google Scholar] [CrossRef] [PubMed]
  6. Wang, L.; Wang, Y.; Ye, D.; Liu, Q. Review of the 2019 novel coronavirus (SARS-CoV-2) based on current evidence. Int. J. Antimicrob. Agents 2020, 55, 105948, Erratum in Int. J. Antimicrob. Agents 2020, 56, 106137. [Google Scholar] [CrossRef] [PubMed]
  7. WHO, World Health Organization. WHO Director-General’s Opening Remarks at the Media Briefing on COVID-19. 11 March 2020. Available online: https://www.who.int/director-general/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19---11-march-2020 (accessed on 29 December 2021).
  8. Cappi, R.; Casini, L.; Tosi, D.; Roccetti, M. Questioning the seasonality of SARS-COV-2: A Fourier spectral analysis. BMJ Open 2022, 12, e061602. [Google Scholar] [CrossRef] [PubMed]
  9. Davis, J.T.; Chinazzi, M.; Perra, N.; Mu, K.; y Piontti, A.P.; Ajelli, M.; Dean, N.E.; Gioannini, C.; Litvinova, M.; Merler, S.; et al. Cryptic transmission of SARS-CoV-2 and the first COVID-19 wave. Nature 2021, 600, 127–132. [Google Scholar] [CrossRef]
  10. Lemos-Paião, A.P.; Silva, C.J.; Torres, D.F.M. A new compartmental epidemiological model for COVID-19 with a case study of Portugal. Ecol. Complex. 2020, 44, 100885. [Google Scholar] [CrossRef]
  11. Ndaïrou, F.; Area, I.; Nieto, J.J.; Silva, C.J.; Torres, D.F.M. Fractional model of COVID-19 applied to Galicia, Spain and Portugal. Chaos Solitons Fractals 2021, 144, 110652. [Google Scholar] [CrossRef]
  12. Silva, C.J.; Cruz, C.; Torres, D.F.M.; Muñuzuri, A.P.; Carballosa, A.; Area, I.; Nieto, J.J.; Fonseca-Pinto, R.; Passadouro, R.; Soares dos Santos, E.; et al. Optimal control of the COVID-19 pandemic: Controlled sanitary deconfinement in Portugal. Sci. Rep. 2021, 11, 3451. [Google Scholar] [CrossRef] [PubMed]
  13. Tang, Y.; Serdan, T.D.A.; Alecrim, A.L.; Souza, D.R.; Nacano, B.R.M.; Silva, F.L.R.; Silva, E.B.; Poma, S.O.; Gennari-Felipe, M.; Iser-Bem, P.N.; et al. A simple mathematical model for the evaluation of the long first wave of the COVID-19 pandemic in Brazil. Sci. Rep. 2021, 11, 16400. [Google Scholar] [CrossRef] [PubMed]
  14. Zine, H.; Boukhouima, A.; Lotfi, E.M.; Mahrouf, M.; Torres, D.F.M.; Yousfi, N. A stochastic time-delayed model for the effectiveness of Moroccan COVID-19 deconfinement strategy. Math. Model. Nat. Phenom. 2020, 15, 50. [Google Scholar] [CrossRef]
  15. Hernandez-Vargas, E.A.; Giordano, G.; Sontag, E.; Chase, G.J.; Chang, H.; Astolfi, A. Second special section on systems and control research efforts against COVID-19 and future pandemics. Annu. Rev. Control 2021, 51, 424–425. [Google Scholar] [CrossRef]
  16. Agarwal, P.; Nieto, J.J.; Ruzhansky, M.; Torres, D.F.M. Analysis of Infectious Disease Problems (COVID-19) and Their Global Impact, Infosys Science Foundation Series in Mathematical Sciences; Springer: Singapore, 2021. [Google Scholar] [CrossRef]
  17. Arino, J. Describing, modelling and forecasting the spatial and temporal spread of COVID-19: A short review. In Mathematics of Public Health; Fields Institute Communications; Springer: Cham, Switzerland, 2022; Volume 85, pp. 25–51. [Google Scholar] [CrossRef]
  18. Wang, L.; Zhang, Q.; Liu, J. On the dynamical model for COVID-19 with vaccination and time-delay effects: A model analysis supported by Yangzhou epidemic in 2021. Appl. Math. Lett. 2022, 125, 107783. [Google Scholar] [CrossRef]
  19. Arino, J.; van den Driessche, P. Time delays in epidemic models, modeling and numerical considerations. Delay Differ. Equ. Appl. 2006, 13, 539–578. [Google Scholar] [CrossRef]
  20. Silva, C.J.; Maurer, H. Optimal control of HIV treatment and immunotherapy combination with state and control delays. Optim. Control Appl. Meth. 2020, 41, 537–554. [Google Scholar] [CrossRef]
  21. Silva, C.J.; Maurer, H.; Torres, D.F.M. Optimal control of a tuberculosis model with state and control delays. Math. Biosci. Eng. 2017, 14, 321–337. [Google Scholar] [CrossRef]
  22. Tipsri, S.; Chinviriyasit, W. The effect of time delay on the dynamics of an SEIR model with nonlinear incidence. Chaos Solitons Fractals 2015, 75, 153–172. [Google Scholar] [CrossRef]
  23. Fine, P.E. The interval between successive cases of an infectious disease. Am. J. Epidemiol. 2003, 158, 1039–1047. [Google Scholar] [CrossRef]
  24. He, X.; Lau, E.H.Y.; Wu, P.; Deng, X.; Wang, J.; Hao, X.; Lau, Y.C.; Wong, J.Y.; Guan, Y.; Tan, X.; et al. Temporal dynamics in viral shedding and transmissibility of COVID-19. Nat. Med. 2020, 26, 672–675. [Google Scholar] [CrossRef]
  25. Xin, H.; Li, Y.; Wu, P.; Li, Z.; Lau, E.H.Y.; Qin, Y.; Wang, L.; Cowling, B.J.; Tsang, T.K.; Li, Z. Estimating the latent period of coronavirus disease 2019 (COVID-19). Clin. Infect. Dis. 2022, 74, 1678–1681. [Google Scholar] [CrossRef]
  26. WHO, World Health Organization. Novel Coronavirus (2019-nCoV): Situation Report-7. Available online: https://apps.who.int/iris/handle/10665/330771 (accessed on 29 December 2021).
  27. Muller, C.P. Do asymptomatic carriers of SARS-COV-2 transmit the virus? Lancet Reg. Health—Europe 2021, 4, 100082. [Google Scholar] [CrossRef]
  28. Peng, L.; Yang, W.; Zhang, D.; Zhuge, C.; Hong, L. Epidemic analysis of COVID-19 in China by dynamical modeling. arXiv 2020, arXiv:2002.06563v2. [Google Scholar] [CrossRef]
  29. Calleri, F.; Nastasi, G.; Romano, V. Continuous-time stochastic processes for the spread of COVID-19 disease simulated via a Monte Carlo approach and comparison with deterministic models. J. Math. Biol. 2021, 83, 34. [Google Scholar] [CrossRef]
  30. Rihan, F.A.; Alsakaji, H.J.; Rajivganthi, C. Stochastic SIRC epidemic model with time-delay for COVID-19. Adv. Differ. Equ. 2020, 2020, 502. [Google Scholar] [CrossRef]
  31. Zaitri, M.A.; Bibi, M.O.; Torres, D.F.M. Optimal control to limit the spread of COVID-19 in Italy. Kuwait J. Sci. 2021, Special Issue on COVID. 1–14. [Google Scholar] [CrossRef]
  32. Lozano, M.A.; Orts, Ò.G.I.; Piñol, E.; Rebollo, M.; Polotskaya, K.; Garcia-March, M.A.; Conejero, J.A.; Escolano, F.; Oliver, N. Open Data Science to Fight COVID-19: Winning the 500k XPRIZE Pandemic Response Challenge. In Machine Learning and Knowledge Discovery in Databases. Applied Data Science Track. ECML PKDD 2021. Lecture Notes in Computer Science; Dong, Y., Kourtellis, N., Hammer, B., Lozano, J.A., Eds.; Springer: Cham, Switzerland, 2021. [Google Scholar] [CrossRef]
  33. Miikkulainen, R.; Francon, O.; Meyerson, E.; Qiu, X.; Sargent, D.; Canzani, E.; Hodjat, B. From Prediction to Prescription: Evolutionary Optimization of Nonpharmaceutical Interventions in the COVID-19 Pandemic. IEEE Trans. Evol. Comput. 2021, 25, 386–401. [Google Scholar] [CrossRef]
  34. Giordano, G.; Blanchini, F.; Bruno, R.; Colaneri, P.; Di Filippo, A.; Di Matteo, A.; Colaneri, M. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat. Med. 2020, 26, 855–860. [Google Scholar] [CrossRef]
  35. Liu, Z.; Magal, P.; Seydi, O.; Webb, G. A COVID-19 epidemic model with latency period. Infect. Dis. Model. 2020, 5, 323–337. [Google Scholar] [CrossRef]
  36. van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef]
  37. Lemos-Paião, A.P.; Silva, C.J.; Torres, D.F.M. A cholera mathematical model with vaccination and the biggest outbreak of world’s history. AIMS Math. 2018, 3, 448–463. [Google Scholar] [CrossRef]
  38. Rogers, J.W. Locations of roots of polynomials. SIAM Rev. 1983, 25, 327–342. [Google Scholar] [CrossRef]
  39. Kuang, Y. Delay Differential Equations with Applications in Population Dynamics; Mathematics in Science and Engineering, 191; Academic Press, Inc.: Boston, MA, USA, 1993. [Google Scholar]
  40. Niculescu, S.-I. Delay Effects on Stability; Lecture Notes in Control and Information Sciences, 269; Springer: London, UK, 2001. [Google Scholar]
  41. Shampine, L.F.; Reichelt, M.W. The MATLAB ODE suite. SIAM J. Sci. Comput. 1997, 18, 1–22. [Google Scholar] [CrossRef]
  42. Silva, C.J.; Cantin, G.; Cruz, C.; Fonseca-Pinto, R.; Passadouro, R.; Soares dos Santos, E.; Torres, D.F.M. Complex network model for COVID-19: Human behavior, pseudo-periodic solutions and multiple epidemic waves. J. Math. Anal. Appl. 2022, 514, 125171. [Google Scholar] [CrossRef] [PubMed]
  43. United Nations. The 2022 Revision of World Population Prospects. 2022. Available online: https://population.un.org/wpp/ (accessed on 29 December 2021).
  44. Cheynet, E. Generalized SEIR Epidemic Model (Fitting and Computation); Transform to Open Science, NASA; Zenodo: Washington, DC, USA, 2020. [Google Scholar] [CrossRef]
Figure 1. Dynamics of model (3) with τ = 3 days. Parameter values as in Table 1.
Figure 1. Dynamics of model (3) with τ = 3 days. Parameter values as in Table 1.
Axioms 11 00400 g001
Figure 2. Dynamics of model (3) with τ [ 0 , 6 ] days. Parameter values as in Table 1.
Figure 2. Dynamics of model (3) with τ [ 0 , 6 ] days. Parameter values as in Table 1.
Axioms 11 00400 g002
Figure 3. Predictions for Italy from the delayed model (24) with τ = 0 and u { 0 % , 20 % , 40 % , 60 % } , between 18 October 2020, and 19 January 2021.
Figure 3. Predictions for Italy from the delayed model (24) with τ = 0 and u { 0 % , 20 % , 40 % , 60 % } , between 18 October 2020, and 19 January 2021.
Axioms 11 00400 g003
Figure 4. Predictions for Italy from the delayed model (24) with u = 20 % and τ { 0 , 3 , 6 } days, between 18 October 2020, and 19 January 2021.
Figure 4. Predictions for Italy from the delayed model (24) with u = 20 % and τ { 0 , 3 , 6 } days, between 18 October 2020, and 19 January 2021.
Axioms 11 00400 g004
Table 1. Parameter values used in the simulations of Section 4.1.
Table 1. Parameter values used in the simulations of Section 4.1.
ParameterValueUnitsRef
b1 Assumed
μ 1 Assumed
δ 1 Assumed
α 1day 1 Assumed
β 12day 1 Assumed
γ 1day 1 Assumed
λ 1day 1 Assumed
t f 30dayAssumed
Table 2. Parameter values used in the simulations of Section 4.2, modeling the spread of the epidemic of COVID-19 in Italy for the period of three months starting 18 October 2020.
Table 2. Parameter values used in the simulations of Section 4.2, modeling the spread of the epidemic of COVID-19 in Italy for the period of three months starting 18 October 2020.
ParameterValueUnitsRef.
b7.391‰ [43]
μ 10.658‰ [43]
α 1.1775day 1 [44]
β 3.97day 1 [44]
γ 0.0048day 1 [44]
λ 0.0182256day 1 [44]
δ 0.1432 [44]
t f 90dayAssumed
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zaitri, M.A.; Silva, C.J.; Torres, D.F.M. Stability Analysis of Delayed COVID-19 Models. Axioms 2022, 11, 400. https://doi.org/10.3390/axioms11080400

AMA Style

Zaitri MA, Silva CJ, Torres DFM. Stability Analysis of Delayed COVID-19 Models. Axioms. 2022; 11(8):400. https://doi.org/10.3390/axioms11080400

Chicago/Turabian Style

Zaitri, Mohamed A., Cristiana J. Silva, and Delfim F. M. Torres. 2022. "Stability Analysis of Delayed COVID-19 Models" Axioms 11, no. 8: 400. https://doi.org/10.3390/axioms11080400

APA Style

Zaitri, M. A., Silva, C. J., & Torres, D. F. M. (2022). Stability Analysis of Delayed COVID-19 Models. Axioms, 11(8), 400. https://doi.org/10.3390/axioms11080400

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop