Stability Analysis of Delayed COVID-19 Models

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.


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.
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 arXiv:2208.06649v1[q-bio.PE] 13 Aug 2022 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 SEIQRP mathematical model for COVID-19.Considering the normalized model of the delayed SEIQRP 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 seiqrpw 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.

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 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: where the state variables are subject to the initial conditions , and P(0) = P 0 .

The Normalized seiqrp 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 Ṅ(t) = (b − µ)N(t).Therefore, the normalized seiqrp delayed model is given by (3) The state variables for system (3) are subject to the following initial conditions: , 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.

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: from which the disease free equilibrium, Σ 0 , is given by while the endemic equilibrium point, Σ + , is given by with Following the method of van den Driessche [36], one easily compute the following basic reproduction number: 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].

Stability of the Normalized seiqrp 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: and x 6 (t) = p(t) − p, where (s, r, ī, q, r, p) denotes any equilibrium point of system (3).The linearized system of (3) takes the form where , and .
The characteristic equation of system (3), for any equilibrium point, is given by We are now in a position to prove our first two results.
Theorem 1 (Stability of the disease free equilibrium of system (3)).If R 0 < 1, then the disease free equilibrium Σ 0 is locally asymptotically stable for any time-delay τ ≥ 0. If R 0 > 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 where We divide the proof into the non-delayed and delayed cases.
(i) Let τ = 0.In this case, the Equation (10) becomes 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, 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 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 By adding up the squares of both equations, and using the fundamental trigonometric formula, we obtain that which is equivalent to 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) (13) 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: and where Proof.The characteristic Equation ( 9), computed at the endemic equilibrium point Σ + , is given by where (i) Let τ = 0.In this case, the Equation ( 16) becomes where ∆1 = L 1 + L1 , ∆2 = L 2 + L2 and ∆3 = L 3 + L3 .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 Using the Routh-Hurwitz criterion [38], we know that all roots of P3 (y, 0) have negative real parts if, and only if, the coefficients of P3 (y, 0) are strictly positive and ∆ (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 By adding up the squares of both equations, and using the fundamental trigonometric formula, we obtain that where Assume that the basic reproduction number R 0 satisfies relations ( 14) and ( 15) with the following condition: Then, In contrast, if R 0 satisfies relations ( 14) and ( 15) with the condition then we have Thus, Under the assumption that the basic reproduction number R 0 satisfies relations ( 14) and ( 15), we have Therefore, if we assume that the basic reproduction number R 0 satisfies relations ( 14) and (15) with condition (20), then if R 0 satisfies relations ( 14) and (15) with condition ( 19), then we have and also equivalent to 1 − (R 0 − 2) 2 > 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).

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 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: where the total population N(t) is given by The state variables are subject to the following initial conditions: = 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.

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 The disease free equilibrium of model ( 24), Σ 1 , is given by while the endemic equilibrium point for system (24), Σ + V , is given by where Following the method from van den Driessche [36], we obtain the following basic reproduction number, denoted by R0 :

Stability of the Normalized seiqrpw Delayed Model with Vaccination
Consider the following coordinate transformation: where (s, r, ī, q, r, p, w) denotes an equilibrium point of system (24).The linearized system of (24) takes the form Ẋ(t) = Ã0 X(t) + Ã1 X(t − τ), (28) where X = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , The characteristic equation of system ( 24) is given by We are also able to prove stability results for the normalized seiqrpw delayed model with vaccination.
Theorem 3 (Stability of the disease free equilibrium of system (24)).If R0 < 1, then the disease free equilibrium Σ 1 is locally asymptotically stable for any time-delay τ ≥ 0. If R0 > 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 where (i) Let τ = 0.In this case, the Equation ( 30) becomes 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 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, Therefore, we have proved that the disease free equilibrium, Σ 1 , is locally asymptotically stable for τ = 0, whenever R0 < 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 By adding up the squares of both equations and using the fundamental trigonometric formula, one has which is equivalent to Hence, we have w 2 < 0, which is a contradiction.Therefore, we have proved that if R0 < 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 R0 > 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 have negative real parts.It is easy to see that q(0) = Γ 2 (0) < 0, because we are assuming R0 > 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 R0 > 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 R0 satisfies the following relations: where Proof.The characteristic Equation ( 29), computed at the endemic equilibrium Σ + V , is given by where (i) Let τ = 0.In this case, Equation (37) becomes where 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 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 R0 > 1, then (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 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 Assume that the basic reproduction number R0 satisfies relations (34) and (35) with the condition Then, In contrast, if R0 satisfies relations (34) and (35) under the condition then we have Thus, K * 1 > 0. Under the assumption that the basic reproduction number R0 satisfies relations (34) and (35), we have Therefore, if we assume that the basic reproduction number R0 satisfies relations (34) and (35) with condition (41), then if R0 satisfies (34) and (35) with condition ( 40), then we have and also equivalent to 1 − R0 − 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 R0 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).

Numerical Simulations and Discussion
In this section we investigate, numerically, the local stability of the normalized seiqrp and seiqrpw models, illustrating our results from Sections 2 and 3.All numerical computations were performed in the numeric computing environment MATLAB R2019b using the medium order method and numerical interpolation [41].
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].
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).

Table 1 :
Parameter values used in the simulations of Section 4.1.

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.