# Analysis and Parametrical Estimation with Real COVID-19 Data of a New Extended SEIR Epidemic Model with Quarantined Individuals

Authors: C. Legarreta; S. Alonso-Quesada; M. De la Sen
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

During the recent COVID-19 pandemic, quarantine and testing policies have been of vital importance since the causative agent has been a novel virus and no vaccine was developed at the time. In this work, a new epidemiological deterministic model is proposed, analyzed, and discussed. Such a model includes quarantine periods of people with symptoms that have been tested positive, and it will trigger a trace of their close contact, who will be also tested and put in quarantine if the result is positive. Moreover, how the model parameters affect its stability is analyzed with the basic reproduction numberR0. Since the COVID-19 outbreak in Spain (approximately 13/03/2020) until 25/04/2021, different restrictions have been applied. For discussion of a real case study, data have been gathered and used from the Spanish Autonomous Community of Cantabria to estimate the parameters and to see how the restrictions have affected their values. In the parameter estimation process, it has been assumed that the constructed model follows the structure of an ARX model. Finally, by considering that the gathered data are subject to certain errors, the paper discusses how to adequate the model usefulness for its use in fitting and processing data through an estimation mechanism involving the provided daily total and positive performed tests.


---

## Body

## 1. Introduction

In late December of 2019, China reported to the World Health Organization (WHO) several cases of pneumonia apparently linked to the Hunan Wholesale Seafood Market in Wuhan City, and they pointed out that this disease might be caused by a new virus. Chinese scientists analyzed the virus from a hospitalized person, and the 11th of January, China confirmed that the abnormal arising of pneumonia cases was due to a novel coronavirus, similar to SARS and MERS, so they named it SARS-CoV-2 (afterwards it will be renamed to COVID-19). It was also found out that fever and sore cough were the common symptoms, and that pneumonia was a less common consequence of the disease causing virus [1]. An abundant background literature on many medical aspects and transmission mechanisms of COVID-19 is available, some of them being linked with general description of epidemic models. See, for instance, [2–23] and references there in and also [24–28] for some general results on epidemic model parameterizations of interest. Some specific related details are now described. By late January of 2020, some cities of China locked down so as to slow down the increase in hospitalized people, and some airlines decided to suspend flights to (or from) China. On the 30th of January, the WHO (World Health Organization) reported that there were 7,818 confirmed cases all over the world, from which 7,736 cases were in China, 1,370 were severe and 170 deaths [21]. Moreover, the Centers for Disease Control and Prevention (CDC) wrote a press release which outlined the first case of person-to-person spread in U.S, and consequently world countries geared up [8]. The applied prevention measures were not enough to avoid the virus spreading and on March 11 determined that COVID-19 was a global pandemic [19] since there were 132,758 confirmed cases in the world and 4,955 deaths [18]. On March 17, Europe and Russia closed the borders, and by the time, many countries announced lock downs (for instance, Spain, France, and so on) [1].When Spain imposed a quarantine period on the 13th of March of 2020, its epidemiological situation was critical: there were 4,209 confirmed cases from which 575 were new ones, and 42% of the total cases were hospitalized [5]. Even if the rules were very restrictive (only on necessary situations it was allowed going outside and wearing the face mask was mandatory), the cases continued raising until 25/03/2020 with a pick value of 720 new cases and 825 deaths [6], and it was not until 25/05/2020, when the number of new cases was approximately 2,000 and 68 the number of deaths [6], that more permission was given to the inhabitants (i.e., it was possible working out outside at certain time ranges). On the 21th of June 2020, the Spanish Government announced the so-called “new normality,” the mask face and social distancing (1.5 m) were still mandatory although at this point each autonomous community could impose more constrains but their capability was limited [7]. During summer vacation, the situation maintained stable; however, when people started working (approximately after 09/09/2020), the epidemiological situation worsen. To tackle this problem, each Autonomous Community (AC) combined and applied differently restrictions (as, for instance, curfew or bars close down), and considering the later results, inhabitants cast doubt on whether the imposed rules were worth or not. Therefore, a mathematical model which could represent the COVID-19 behavior over certain conditions could be very useful.The aim of epidemiological models is to understand and predict the dynamics of the propagation of a disease, and it must exhibit an equilibrium between three elements [11], namely, accuracy, transparency, and flexibility. Accuracy is the capacity to reproduce the real life behavior and a necessary feature to establish good control actions, perhaps the complexity of the model is proportional to the mentioned accuracy. On the other hand, transparency considers how the model components interact in the dynamics, and it has to be noticed that it is a property that is opposed to complexity and consequently to accuracy. Finally, flexibility is the ability to adapt the model to other circumstances.If the model fits well with the real disease propagation dynamics, it can be used to infer optimal control actions, so it makes possible implementing more efficient measures and optimizing the use of resources. In 2001, the epidemic foot-and-mouth outbreak in the UK and the result of three different models applied to this case predicted a large-scale spreading and a better control of the epidemic if the infected and exposed animals were culled [26]. In many research studies, it has been considered a spreading of a very high death rate virus such as smallpox [11] to verify whether a massive vaccination, which could generate health problems, is advantageous or not. However, the obtained results differ from each other since there are many epidemiological uncertainties. As far as COVID-19 is a pandemic, there are many data sources that provide information about the epidemiological situation, so it is feasible to create a parametric model from real data and reduce epidemiological uncertainties. Thus, clear control strategies could be built up and implemented.In this work, a deterministic mathematical model is built up, which spans the well-known SEIR epidemic model by adding people in quarantine and those being asymptomatic. To reduce the model dimension, the infected people with symptoms and asymptomatic will be gathered, so a SEIQR epidemic model is obtained which is integrated by the following subpopulations: Susceptible (S), Exposed (E), Infectious (I), Quarantined (Q), and Recovered or Immune (R). In this model, it is assumed that only exposed and infectious people will cause new infections since people in Q, which represent the individuals who stay in hospital and at home during some period due to the disease, are concerned about their illness, and they try not being in contact with others. On the other hand, it will be considered that some percentage of the individuals in the infectious subpopulation I, more precisely the symptomatic ones, will be reported and carried to Q, and it will trigger a tracing method as it is shown in [3, 23]. Thus, we can say that there are the subpopulations including infected individuals which are the exposed, the infectious, and the quarantined subpopulations. In particular, the asymptomatic infected subpopulation is included in the exposed one.After constructing the SEIQR model schematic, defining how individuals interact, and defining the necessary parameters, the respective nonlinear system of differential equations will be obtained and analyzed. The theoretical expression for the basic reproduction number, which is a feature dependent on the system’s parameters and a threshold that determines whether the disease is out of control or not, will be obtained through the next generation matrix [27]. For stability analysis purposes, the equilibrium points, that is, the disease-free one and the endemic equilibrium one, which determine the states in which the disease is extinguished from the society and when the disease prevail (i.e., the usual flu disease), respectively, are calculated. The system is linearized around the equilibrium points, and the eigenvalue problem is solved to evaluate when their real parts become positive/negative (condition instability/stability of the system).Some of the model parameters (for instance, the virus incubation period and the probability to get infected by symptomatic/asymptomatic people) are a disease characteristic, so their values have been collected via other studies [2, 3, 12, 13, 27, 29]. On the other hand, in [30], a model with susceptible, exposed, infectious, quarantined, and removed is discussed with parametrical values fitted to Saudi Arabia data which test the model and analyze the effects of hospital beds availability, quarantines, and media effects on the foreseen dynamics. In [31], a sliding-mode vaccination controller is designed to fight against COVID-19 and testing simulations are performed on data of Iran and Russia. In [32], machine learning designs are used in a susceptible-infected-recovered model to discuss the influence of vaccination in the COVID-19 pandemic evolution in Brazil as well as the impacts of the immunization speed and vaccines efficacy.However, parameters that depend on the society behavior (i.e., the average number of contacts per day) must be determined with the respective real data. Since the autonomous community of Cantabria gathers wide data about the epidemiological situation, it has been used to estimate the unknown parameters, such as the average contact rate and the tracing effectiveness. It has been assumed that the discrete differential expressions of variables such as Q and I follow an ARX model structure. Thus, the ARX model parameters have been estimated with the least square method and their values have been used to infer the unknown epidemiological values of the parameters. Finally, the differential equations have been solved with the obtained parameters, and the results have been compared with real data.The paper is organized as follows. Section2 describes the new proposed epidemic SEQIR (Susceptible-Exposed-Quarantined-Infectious-Recovered) model. Such a model has several specific novel characteristics, whose effects are integrated and analyzed together in a coordinated way, namely,(a)
The proportions of asymptomatic and symptomatic individuals in the quarantined group and in the infectious group have, in general, different average transmission rates of contagion to susceptible individuals as it is also potentially distinct from the above mentioned ones the average transmission rate of the exposed individuals to the susceptible ones.(b)
Various average transmission rates are parameterized by the infective contacts between individuals. In particular, the values of the average transmission rates are directly proportional to the average potentially infectious contacts. Those contacts occur among the various infected types of individuals (that is, those being currently allocated in the exposed subpopulation, including the asymptomatic infective, in the symptomatic infectious subpopulation, and in the quarantined subpopulation) with susceptible individuals.(c)
Various average transmission rates are also parameterized in a direct proportion by the corresponding probabilities of infecting susceptible individuals.(d)
The model is also parameterized by the average rate of attendance to the doctor for disease positivity testing or evaluation of symptoms checking and by the effectiveness in tracing the positive cases.(e)
The increase in the attendance rate of symptomatic suspect individuals contributes positively to the incorporation to the quarantined subpopulation of those being positively tested.(f)
The total day-to-day quarantined subpopulation incorporates the individuals who have been tested positive. This includes the hospitalized individuals including those under intensive care and also those quarantined at their homes having no serious symptoms.The positivity of the model, in terms of the nonnegativity of all the components of the state trajectory solution under any given finite nonnegative initial conditions, is also proved and discussed in that section. Section2 also includes a brief discussion on other approaches given for the pandemic study including the potential vaccination and treatment control relevance in the damping of the infection as well as the influences of the intervention measures in the pandemic evolution. On the other hand, Section 3 deals with the disease-free and the endemic equilibrium points which are proved to be unique. The disease-free equilibrium point is proved to be locally asymptotically stable if the basic reproduction number is less than unity, while the endemic one is locally asymptotically stable if such a value exceeds unity. In addition, the endemic equilibrium point is proved to be nonreachable in the first case since its infective components are negative so that the eventual convergence of the state trajectory solution to it would contradict the previously proved positivity of the model performed in the former section. Furthermore, it is proved that only one of the two equilibrium points is a global asymptotic attractor depending again on the value of the basic reproduction number. Section 4 gives some simulations with the proposed model being evaluated for the COVID-19 pandemic with demographic data from Spain. Section 5 develops in more detail a case study against COVID-19 which is performed with official data acquisition taken from the Cantabria Autonomous Community (a region in North Spain). In particular, the estimations of the average periods lasted by people to get to a Health Center or to a doctor for eventual positivity inspection, the efficiency of the positivity tracing, and the average number of contagion contacts are estimated from real data along the initial confinement period and also for a subsequent period. In this subsequent period, the formerly applied strong public intervention measures, of almost total general confinement of the whole population, were left while still public intervention rules were maintained, for instance, a duty of general use of masks, social distance keeping, partial restrictions of mobility, closing of nocturnal amusement places or very significant limitation rules on the number of attendees to spectacles and restaurants, and so on. Finally, conclusions end the paper. A subsection with a glossary of the main symbols is listed below, and the proofs of the stated mathematical results are given in appendixes.


### 1.1. Nomenclature

S, E, Q, I, R: subpopulations of susceptible, exposed, quarantined, infectious, and recovered individuals. It is assumed that the exposed subpopulation E includes the asymptomatic infective individuals, while the symptomatic infective (or infectious) integrated in the subpopulation I.N: total populationΛ: average recruitment parameterβe, βa, βs: exposed to susceptible, asymptomatic to susceptible, and symptomatic to susceptible average transmission ratesσ: average incubation rateγ: average recovery rateρ: average immunity loss rateδ: average quarantine rateμ: average natural death rateα: average disease mortality ratep1: average probability of evolution from exposed to symptomaticp2: average probability of evolution from exposed to asymptomaticλ: average rate of symptomatic visits to doctor for testingce=βe/pe, ca=βa/pa, cs=βs/ps: average numbers of close contacts per day of a member from a group of (either exposed, asymptomatic or symptomatic infectious) to susceptible; and pe, pa, and ps are the respective probabilities to infect a susceptible individualt,p: daily total and positive testsr: average number of recovered quarantined individuals who were confined because positive testing


## 1.1. Nomenclature

S, E, Q, I, R: subpopulations of susceptible, exposed, quarantined, infectious, and recovered individuals. It is assumed that the exposed subpopulation E includes the asymptomatic infective individuals, while the symptomatic infective (or infectious) integrated in the subpopulation I.N: total populationΛ: average recruitment parameterβe, βa, βs: exposed to susceptible, asymptomatic to susceptible, and symptomatic to susceptible average transmission ratesσ: average incubation rateγ: average recovery rateρ: average immunity loss rateδ: average quarantine rateμ: average natural death rateα: average disease mortality ratep1: average probability of evolution from exposed to symptomaticp2: average probability of evolution from exposed to asymptomaticλ: average rate of symptomatic visits to doctor for testingce=βe/pe, ca=βa/pa, cs=βs/ps: average numbers of close contacts per day of a member from a group of (either exposed, asymptomatic or symptomatic infectious) to susceptible; and pe, pa, and ps are the respective probabilities to infect a susceptible individualt,p: daily total and positive testsr: average number of recovered quarantined individuals who were confined because positive testing


## 2. Model Setting and Its Description and Positivity

Epidemic models can operate in either deterministic or in stochastic frameworks depending on how the population of one block transforms into the others; that is, deterministic models consider that the rates (namely, the chance to be infected) maintain constant with respect to time, and therefore there is a unique solution for each initial condition. In contrast, stochastic models add random or probabilistic variables into the rates, so the model will give a set of probable solutions. Even if human behavior and infections have stochastic components, when the real system (stochastic) is composed by a large group of people, deterministic models could be used to represent the system. Then, the ratios are constant and could be determined by statistical results [11].The graphical representation of the system is made by a group of blocks (or compartments), each one represents a certain part of the population, and the arrows make the connection between two blocks (A ⟶ B) and specify the way in which part of the population from A transforms into the inhabitants of the block B. To build up the model flow chart, it is important to know the features of the tackling disease which are the stages of the disease, incubation period, immunity period, and so on.Regarding COVID-19, the WHO reported that there are two main ways in which the disease is transmitted: presymptomatic and asymptomatic transmissions [20]. Presymptomatic transmission, which is possible during the incubation period (time between exposure and the symptoms onset) which is 5-6 days [17, 20], is more likely to happen 1–3 days before symptoms appear [26]. The preliminary way of transmission is via symptomatic cases, but even if the percentage of asymptomatic cases is 16%-17%, [2, 13], it has been seen that they show similar viral loads [10, 13, 26, 27, 29]. Some research studies have spotlighted the importance of presymptomatic and asymptomatic cases since they contribute to the virus high spread [27, 29]. Therefore, many countries were forced to include quarantine periods to reduce the spreading. The flow chart depicted in Figure 1 shows a deterministic model of COVID-19 which includes the characteristics mentioned before. Note that Figure 1(a) shows five different compartments, namely,(i)
Susceptible (S): this class is constructed by a group of individuals who can get the disease and are not yet infected.(ii)
Exposed (E): those individuals of a group S who have been in contact with infected people and therefore are infected too, but they do not present symptoms because they are still incubating the virus. They move freely; thus, they will maintain an average number of close contacts per day cfree.(iii)
Symptomatic (Sy): this group consists of people who got the disease and present symptoms. Soon after the incubation time, an individual of Sy will probably move freely (the average number of close contacts per day will be equal to cfree) since the host manifests presymptoms, which commonly are confused with tiredness or symptoms of a common cold/flu. After the list of symptoms will become larger and their intensity will increase, therefore, he/she will go to the doctor and put in quarantine. It is assumed that the average number of close contacts of people from Sy is lower than cfree.(iv)
Quarantined (Q): this group is integrated by subgroups coming from S and Sy. People from Sy who presented clear symptoms will have tested positive and put in quarantine. Their positive results will trigger a tracing through their close contacts, which is assumed that are individuals from S. Traced contacts will be tested, and those with positive result will be put in quarantine. It is assumed that the average number of contacts per day of people in Q is nearly zero because they are aware of their disease and they will interact less with others. Note that Q covers people who have tested positive and are staying at home or hospital; therefore, only people from Q might die due to the disease.(v)
Asymptomatic (A): a group of infected people who does not present symptoms. They move freely with an average number of close contacts equal to cfree.(vi)
Recovered (R): when infected, people recover from the disease and have some immunity to it during a period of time.Figure 1
(a) Flow chart of the epidemic model and (b) reduced flow chart of the epidemic model.
(a)(b)The newborns are introduced in the classS with a recruitment Λ, and from the classes, it removed a portion with a rate μ (daily natural death rate). It is assumed that the population is homogeneously mixed, so it is possible to determine the transmission rates βe, βa, and βs, which indicate the probability to infect susceptible people when they are in contact with exposed, asymptomatic, and symptomatic individuals, respectively. The corresponding expressions are βe=cepe, βa=capa, and βs=csps where ci is the average number of close contacts per day of a member from exposed, asymptomatic, and symptomatic, respectively, to susceptible, while pe, pa, and ps are the respective probabilities to infect a susceptible individual. As it was mentioned before, it will be assumed that ce=ca=cfree, ce>cs>0, and pa≈ps>pe. These are necessary parameters to specify how susceptible people get infected; one individual from S has a probability E/N, I/N, and A/N of contacting a person from E, I, and A, respectively, so the chance to get infected is βeE+βsSy+βaA/N. Therefore, all susceptible people who will get disease per day is(1)βeE+βsSy+βaASN.Exposed people after the incubation time (σ−1 days) will evolve to asymptomatic people with a probability p2. Those ones evolving to symptomatic individuals may not consider, at the early beginning, that they are still sick. However, after λ−1 days (λ indicates the rate at which people from group Sy visit the doctor, and it is supposed that its value will be larger as the test stock, sanitary resources, and the disease impact on society are increased), their symptoms will increase and they will go to the doctor, test positive, and put in quarantine. This will start a tracing, which effectiveness is determined by p1, of the positive tested symptomatic people with close contacts csλp1SyS/N, and those who test positive βsλp1SyS/N will go to Q. Therefore, the number of new exposed per day will be equal to the susceptible people who get the disease per day (see expression (1) minus traced people who have tested positive):(2)βeE+βs1−p1λSy+βaASN.People inQ include those who have been tested positive, which means hospitalized ones, people in Intensive Care Unit (ICU), and those who are in their homes. Some people from Q will remain there during the quarantine period (δ−1 days), and others will die from COVID-19 with a rate α. After γ−1 days, both asymptomatic and symptomatic individuals will recover from the disease, and they will be part of the recovered population R during a period ρ−1. Now, six blocks are shown in the whole scheme from Figure 1(a), so that six differential equations are needed to represent the system. Managing with six equations might be a difficult task; therefore, a new compartment I involving together A and Sy is defined to reduce the number of equations from six to five, that is,(3)I=A+Sy,A=p2I,Sy=1−p2I,and in Figure 1(b), the block I (infectious) substitutes the blocks A and Sy. In the following, the dependence on time is omitted in the equations for the sake of exposition and simplicity when no confusion is expected. The differential system representing the deterministic model is given by the following set of first-order differential equations:(4)S˙=Λ+ρR−βeE+p2βa+1−p2βsISN−μS,(5)E˙=βeE+p2βa+1−p21−λp1βsISN−σ+μE,(6)Q˙=1−p2βsp1λISN+λ1−p2I−μ+α+δQ,(7)I˙=σE−λ1−p2+μ+γI,(8)R˙=γI+δQ−ρ+μR,subject to nonnegative initial conditions S0=S0,E0=E0,Q0=Q0,I0=I0, and R0=R0. The above model is being considered in the sequel for analysis. All the parameters are nonnegative and, furthermore, p1,p2∈0,1. The total population and its time derivative become(9)N=S+Q+E+I+R;N˙=Λ−μN−αQ.The following result holds on nonnegativity of the trajectory solution for any given finite nonnegative initial conditions. The proof is given in Appendix A.Theorem 1.
Assume thatx0=x0=S0,E0,Q0,I0,R0T∈R5+. Then, the following properties hold:(i)
xt=St,Et,Qt,It,RtT∈R5+ for all t∈R+ so that the five subpopulations are nonnegative for all time(ii)
The total population is nonnegative and bounded for all time and then the subpopulations are also bounded for all time.Close techniques for evaluating the nonnegativity of the solution have been used in some background literatures for other epidemic models. See, for instance, [27–32]. Also, it can be pointed out that other related alternative mathematical models for infective disease have been proposed in the recent literature with applications to parameterizations related to COVID-19. For instance, in [33], mixed contagions of susceptible individuals from asymptomatic and symptomatic infectious ones are studied based on provided official data taken from some European countries. In [34], SARS-CoV-2 (COVID-19) transmission data from Spain and Italy were analyzed and processed to estimate the transmission rate levels through SIR models with undelayed and delayed propagation infection related to the periods of different public intervention measures. In [35], confinements and quarantines have been considered as impulsive actions which decrease at certain time instants, and along certain duration time periods, the amount of individuals is able to produce contagions. In [36], an analysis of the contagion propagation levels from small amounts of exposed and infected outsiders coming in a certain habitat was studied under eventual delayed resusceptibility. In [37], a model for COVID-19 including a potential herd immunity has been developed for Austria, Luxembourg, and Sweden. Also, the effect of delaying the vaccination while assuming unlimited vaccine units supply is studied in [38]. A very complete recent work on epidemic modeling is given in [39]. Such a well-organized work discusses different deterministic, stochastic, and optimization models and includes also a discussion on the pandemic evolution in France. On the other hand, some mathematical techniques of interest for its application in mathematical epidemic models related to positivity and stability as well as numerical methods, estimation methods, or probability and statistics tools can be found in [39–47].


## 3. Equilibrium Points and Stability Analysis

Now, the disease-free and the endemic equilibrium points of the above model are seen to be unique, and their components are characterized. It is also seen that the endemic equilibrium is not reachable (or, roughly speaking, it does not exist in view of Theorem1 (i)) for small values of the transmission rate since the infective components are negative.


### 3.1. Equilibrium Points

The equilibrium points are found by zeroing the time derivatives of the subpopulations in (4)–(8). The disease-free equilibrium point satisfies Edfe=Qdfe=Idfde=0. Looking at (8), it is seen that those constraints imply that Rdfe=0 and, from the first one, Sdfe=Ndfe=Λ/μ so that xdfe=Λ/μ,0,0,0,0T is indeed unique. The interpretation is that, in the absence of disease, the whole population is susceptible. Note that, in the absence of recruitment, i.e., if Λ=0, then the disease-free equilibrium point implies the population extinction.The existence, uniqueness, and characterization of the endemic equilibrium point are discussed in the two following results which are, respectively, proved in Appendixes B and C.Lemma 1.
Assume that(10)βc=σ+μλ1−p2+μ+γλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.Then, the following properties hold:(i)
We have(11)Iee=fIEee,fI=σλ1−p2+μ+γ.(ii)
We have(12)See=gSNee,gS=σ+μλ1−p2+μ+γλ1−p2+μ+γβe+p2βa+1−p21−p1λβsσ.Let us define the relative transmission rates related to some transmission rate reference valueβ leading to βe=βerβ, βa=βarβ, and βs=βsrβ. Then, a necessary condition for the endemic equilibrium point to exist is that β>βc.(iii)
The quarantined value at the endemic equilibrium pointQee, if it exists (i.e., if Qee>0), satisfies the subsequent constraint:(13)λ1−p2Ieeμ+α+δ<μ1−p2βsp1λSeeμ+α+δΛ−αQee+λ1−p2μ+α+δIee=Qee,andQee>0 is unique if λ∈0,λc for some λc≤1/p1. Furthermore, Qee<Λ/μ+α and Iee<Λμ+α+δ/λ1−p2μ+α.(iv)
Qee=fQEee, where(14)fQ=1−p2λσμ+α+δλ1−p2+μ+γ1+βsrp1σ+μλ1−p2+μ+γλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.(v)
We have(15)Ree=γIee+δQeeρ+μ=fREee,where(16)fR=γfI+δfQρ+μ=1ρ+μγfI+1−p2λδμ+α+δfI+βsrp1σ+μσλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.Note that the reference transmission rateβ is fixed for each given triple of transmission rates βe, βa, and βs with their relative values being βer=βe/β, βar=βa/β, and βsr=βs/β.The following brief discussion concludes that a unique feasible endemic equilibrium point can also exist irrespective of the parameterλ which describes the attendance rate to the doctor or Health Service for checking or treatment.Lemma1 leads directly to the next result whose proof is given in Appendix C.Theorem 2.
Assume thatλ∈0,λc and β>βc. Then, the endemic equilibrium point is unique subject to an endemic equilibrium population given by(17)Eee=1−gSΛμ+fQμ+α1−gS+μfI+fR.And, as a result, there is a unique reachable endemic equilibrium point given by (11), (12), (14), (15), and (17). If β=βc, then the endemic equilibrium point coincides with the disease-free one.Remark 1.
Theorem2 establishes that the endemic equilibrium point does not exist for sufficiently small transmission rates so that, in this case, the disease-free equilibrium point is the only feasible equilibrium point. For instance, note that if the three transmission rates are zero, then the susceptible individuals at the endemic equilibrium would be infinity what contradicts Theorem 1(ii) so that, for zero transmission rates, the endemic equilibrium does not exist, as expected.Remark 2.
Note that it becomes attractive to refer the three defined transmission rates to a prescribed reference oneβ by using the relative values of the defined one with respect to such a reference transmission rate. In particular, it is easy in that way to characterize the existence of the endemic equilibrium point referring to such a reference transmission rate exceeding a threshold which depends on the remaining parameters including the three relative transmission rates. It will be seen later on in Section 3.2 that the necessary condition β>βc for the reachability of the endemic equilibrium point (equivalent to its existence from its necessary positivity for that concern) is also equivalent to the basic reproduction number to fulfill R0>1, which implies in turn the instability of the disease-free equilibrium point. Thus, one will conclude that if the disease-free equilibrium point is locally asymptotically stable, then the endemic equilibrium one is not reachable. Simulations corroborate that if R0=1, then the disease-free equilibrium point is locally asymptotically stable. As a result, the disease-free equilibrium point is a global attractor if R0≤1 and equivalently if β≤βc.Remark 3.
Note that the sufficient conditionsλ∈0,λc and β>βc of Theorem 2 for the uniqueness of the endemic equilibrium point are independent of the equilibrium subpopulations amount since it is based on the fact that Iee<Nee.
For the particular case,λ=0 which reflects that people, in average, do not go to the doctor under symptoms since the delayed time to go to the doctor is λ−1=+∞. For such a case, it can be deduced as well that the endemic equilibrium is unique. However, the reasoning has to be specific for the discussion of the particular model (4)–(8) obtained by zeroing the terms affected by a factor λ in the differential equations. Since λ=0 yields to division by zero in some relevant equations in the proof of Lemma 1, the next related result is stated and proved separately from Lemma 1 and Theorem 2. Its proof is given in Appendix D.Theorem 3.
Ifλ=0, then there are a unique disease-free equilibrium point and a unique endemic equilibrium. The first one is xdfe=Λ/μ,0,0,0,0T; similarly for λ≠0, and the endemic one satisfies the subsequent relations:(18)Qee=0;Iee=σEeeμ+γ;Ree=γIeeρ+μ=γσρ+μμ+γEee,(19)See=SeeNeeNee=σ+μμ+γΛμβeμ+γ+p2βa+1−p2βsσ,(20)Eee=Λρ+μμ+γZ−σ+μμ+γZμσ+μρ+μ+γσ+μ+ρ,where(21)Z=βeμ+γ+p2βa+1−p2βsσ.
Furthermore, the endemic equilibrium point is reachable if the reference transmission rate exceeds a minimum threshold according to(22)β>βc0=βcλ=0=σ+μ+μ+γμ+γβer+p2βar+1−p2βsrσ.


### 3.2. Local Asymptotic Stability Analysis

The next technical result, proved in Appendix E, relies on the evaluation of the Jacobian matrix with respect to any equilibrium point. It will be then used to characterize the local asymptotic stability of both equilibrium points.Lemma 2.
The Jacobian matrixJ∗=J∗x∗=∂x˜˙t/∂x˜Tt|x∗ of the linearized trajectory x˜t=xt−x∗ around any equilibrium point x∗=S∗,E∗,Q∗,I∗,R∗T is given by J∗=Jij∗∈R5×5 with entries as follows:(23)J11∗=−μ−βeE∗+p2βa+1−p2βsI∗N∗−S∗N∗2,J12∗=−βeS∗N∗−E∗−p2βa+1−p2βsI∗S∗N∗2,J13∗=βeE∗+p2βa+1−p2βsI∗S∗N∗2,J14∗=−p2βa+1−p2βsS∗N∗−I∗−βeE∗S∗N∗2,J15∗=ρ+βeE∗+p2βa+1−p2βsI∗S∗N∗2,J21∗=βeE∗+p2βa+1−p21−λp1βsI∗N∗−S∗N∗2,J22∗=−σ+μ+βeS∗N∗−E∗−p2βa+1−p21−λp1βsI∗S∗N∗2,J23∗=−βeE∗+p2βa+1−p21−λp1βsI∗S∗N∗2,J24∗=p2βa+1−p21−λp1βsS∗N∗−I∗−βeE∗S∗N∗2,J25∗=−βeE∗+p2βa+1−p21−λp1βsI∗S∗N∗2,J31∗=1−p2βsp1λI∗N∗−S∗N∗2,J32∗=−1−p2βsp1λI∗S∗N∗2,J33∗=−μ+α+δ−1−p2βsp1λI∗S∗N∗2,J34∗=λ1−p2+1−p2βsp1λS∗N∗−I∗N∗2,J35∗=−1−p2βsp1λI∗S∗N∗2,J41∗=J43∗=J45∗=0,J42∗=σ,J44∗=−λ1−p2+μ+γ,J51∗=J52∗=0,J53∗=δ,J54∗=γ,J55∗=−ρ+μ.
The above Jacobian matrix can be decomposed asJ∗=T∗−V with the transmission and transition matrices, respectively, −V (being independent of the transmission rates and of the considered equilibrium point) and T∗=J∗+V, where(24)−V=J∗β=0=−μ000ρ0−σ+μ00000−μ+α+δλ1−p200σ0−λ1−p2+μ+γ000δγ−ρ+μ,by using the relative transmission rates βer=βe/β, βsr=βs/β, and βar=βa/β with respect to the reference transmission rate β.The above result is useful for the local stability analysis around the disease-free and endemic equilibrium points, respectively,xdf and xee, via the respective Jacobian matrices Jdf=J∗xdf=T∗xdf−V and Jee=J∗xdf=T∗xee−V. In the calculation of the entries of the Jacobian matrix around x∗, it has been taken into account that the total population is not constant with respect to time due to the mortality associated to the disease. For instance,(25)∂S/N∂Sx=x∗=∂S/S+E+Q+I+R∂Sx=x∗=N∗×1−S∗×1N∗2=1N∗−S∗N∗2,generates one of the additive terms in J11∗ and so on for the remaining entries of the Jacobain matrix around the considered equilibrium point.The following remark gives some considerations of conceptual interest related to the Jacobain matrix around the equilibrium points and its partition into the sum of the transition and transmission matrices.Remark 4.
Note that the transmission matrix−V is a Metzler stability matrix, i.e., its diagonal entries are negative and its off-diagonal entries are nonnegative so that all its eigenvalues are real negative and its associated fundamental (or transition) matrix e−Vt is nonsingular and has nonnegative entries, for all time. Note also that the transmission matrix T∗ is zero if the transmission rates are zero (i.e., it is zero in the absence of disease). The fact that e−Vt has nonnegative entries for all time implies that any solution trajectory is nonnegative for all time under any given nonnegative initial conditions, in the absence of disease which is a necessary condition for the nonnegativity of the whole model to have a nonnegative solution trajectory under nonnegative initial conditions (Theorem 1).
On the other hand,J∗=−VI−V−1T∗=−I−T∗V−1V is stable since −V is stable if and only if the spectral radius of V−1T∗, identical to that of T∗V−1 is small enough. According to equation (23), this fact happens for the disease-free equilibrium point if β is small enough. It can be proved that this is equivalent also to β<βc and it will be seen later on that this is also further equivalent to the basic reproduction number to be less than one. Finally, according to (10), this also implies that the endemic equilibrium point is not reachable.
Local asymptotic stability and instability results around the equilibrium points are stated in the subsequent result, whose proof is given in Appendix F.Theorem 4.
The following properties hold:(i)
The Jacobian matrix around the disease-free equilibrium point is stable if and only if the basic reproduction number(26)R0=ββc−1=βλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσσ+μλ1−p2+μ+γ<1,with the average critical reference transmission rateβc being defined in (10) and then the disease-free equilibrium point is locally asymptotically stable. Moreover, the spectral radius of V−1Tdf equalizes the basic reproduction number where, such that for the disease-free equilibrium point, the transmission and transition matrices are as follows:(27)Tdf=Tdfβ=Jdf+V=β0−βer0−p2βar+1−p2βsr00βer0λ1−p2p1βsr0000λ1−p21+βsrp100000000000,(28)−V=J∗β=0=−μ000ρ0−σ+μ00000−μ+α+δλ1−p200σ0−λ1−p2+μ+γ000δγ−ρ+μ.(ii)
The disease- free equilibrium point is unstable ifβ>βc⇔R0>1.(iii)
Assume thatλ∈0,λc. Then, the endemic equilibrium point is locally asymptotically stable if and only if β>βc⇔R0>1.(iv)
The above conditions in Properties (ii) and (iii) can be weakened in the sense that the disease-free equilibrium point is locally asymptotically stable if and only ifR0β≤1β≤βc and the endemic one is locally asymptotically stable if and only if R0β≥1β≥βc.Remark 5.
It is now discussed how the calculation of the basic reproduction number can be simplified by manipulating the infective variables, that is, the exposed, infectious, and quarantine subpopulations, only. Note from equations (27)–(28) that the Jacobian matrix at the disease-free equilibrium point is as follows:(29)Jdf=−μ−βe0−p2βa+1−p2βsρ0−σ+μ+βe0p2βa+1−p21−p1λβs000−μ+α+δλ1−p21+βsp100σ0−λ1−p2+μ+γ000δγ−ρ+μ,which makes obvious by its direct inspection that the linearized subsystem of the infective components E, I, and Q around the disease-free equilibrium points is not coupled to the linearized subsystem of noninfective components S and R. Note that the Jacobian matrix of the infective linearized subsystem around the disease-free equilibrium point is given by(30)Jdfi=−VdfiI−Vdfi−1Tdfi=βe−μ+σ0p2βa+1−p21−p1λβs0−μ+α+δλ1−p21+βsp1σ0−λ1−p2+μ+γ,where the corresponding transition and transition matrices are as follows:(31)−Vdfi=−μ+σ000−μ+α+δλ1−p2σ0−λ1−p2+μ+γ;Tdfi=ββer0p2βar+1−p21−p1λβsr00λ1−p2βsrp1000,which concludes that the local stability around the disease-free equilibrium point holds if R0=ρTdfiVdfi<1 which guarantees that Jdfi is a stability matrix since −Vdfi is a stability matrix of eigenvalues −μ+σ, −μ+α+δ, and −λ1−p2+μ+γ. Therefore, the local stability around the disease-free equilibrium point can be analyzed more easily from the Jacobian matrix of the infective linearized subsystem since if Et,It,Qt⟶0 as t⟶∞ for any initial conditions sufficiently close to the disease-free equilibrium point, then the disease-free equilibrium point is locally asymptotically stable and conversely this happens if the eigenvalues of Jdfi have negative real parts.


### 3.3. Global Asymptotic Stability Analysis

The subsequent result extends Theorem4 by proving that the given asymptotic stability results are also global and that only one of the two equilibrium points is a global asymptotically stable attractor depending on the range of values of the basic reproduction number, equivalently, on that of the reference transmission rate for any given fixed triple of relative transmission rates. Its proof is given in Appendix G.Theorem 5.
Assume thatλ∈0,λc. Then, there is only a globally asymptotically stable attractor which is the disease-free equilibrium point if R0β≤1β≤βc and the endemic one if R0>1β>βc.


## 3.1. Equilibrium Points

The equilibrium points are found by zeroing the time derivatives of the subpopulations in (4)–(8). The disease-free equilibrium point satisfies Edfe=Qdfe=Idfde=0. Looking at (8), it is seen that those constraints imply that Rdfe=0 and, from the first one, Sdfe=Ndfe=Λ/μ so that xdfe=Λ/μ,0,0,0,0T is indeed unique. The interpretation is that, in the absence of disease, the whole population is susceptible. Note that, in the absence of recruitment, i.e., if Λ=0, then the disease-free equilibrium point implies the population extinction.The existence, uniqueness, and characterization of the endemic equilibrium point are discussed in the two following results which are, respectively, proved in Appendixes B and C.Lemma 1.
Assume that(10)βc=σ+μλ1−p2+μ+γλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.Then, the following properties hold:(i)
We have(11)Iee=fIEee,fI=σλ1−p2+μ+γ.(ii)
We have(12)See=gSNee,gS=σ+μλ1−p2+μ+γλ1−p2+μ+γβe+p2βa+1−p21−p1λβsσ.Let us define the relative transmission rates related to some transmission rate reference valueβ leading to βe=βerβ, βa=βarβ, and βs=βsrβ. Then, a necessary condition for the endemic equilibrium point to exist is that β>βc.(iii)
The quarantined value at the endemic equilibrium pointQee, if it exists (i.e., if Qee>0), satisfies the subsequent constraint:(13)λ1−p2Ieeμ+α+δ<μ1−p2βsp1λSeeμ+α+δΛ−αQee+λ1−p2μ+α+δIee=Qee,andQee>0 is unique if λ∈0,λc for some λc≤1/p1. Furthermore, Qee<Λ/μ+α and Iee<Λμ+α+δ/λ1−p2μ+α.(iv)
Qee=fQEee, where(14)fQ=1−p2λσμ+α+δλ1−p2+μ+γ1+βsrp1σ+μλ1−p2+μ+γλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.(v)
We have(15)Ree=γIee+δQeeρ+μ=fREee,where(16)fR=γfI+δfQρ+μ=1ρ+μγfI+1−p2λδμ+α+δfI+βsrp1σ+μσλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσ.Note that the reference transmission rateβ is fixed for each given triple of transmission rates βe, βa, and βs with their relative values being βer=βe/β, βar=βa/β, and βsr=βs/β.The following brief discussion concludes that a unique feasible endemic equilibrium point can also exist irrespective of the parameterλ which describes the attendance rate to the doctor or Health Service for checking or treatment.Lemma1 leads directly to the next result whose proof is given in Appendix C.Theorem 2.
Assume thatλ∈0,λc and β>βc. Then, the endemic equilibrium point is unique subject to an endemic equilibrium population given by(17)Eee=1−gSΛμ+fQμ+α1−gS+μfI+fR.And, as a result, there is a unique reachable endemic equilibrium point given by (11), (12), (14), (15), and (17). If β=βc, then the endemic equilibrium point coincides with the disease-free one.Remark 1.
Theorem2 establishes that the endemic equilibrium point does not exist for sufficiently small transmission rates so that, in this case, the disease-free equilibrium point is the only feasible equilibrium point. For instance, note that if the three transmission rates are zero, then the susceptible individuals at the endemic equilibrium would be infinity what contradicts Theorem 1(ii) so that, for zero transmission rates, the endemic equilibrium does not exist, as expected.Remark 2.
Note that it becomes attractive to refer the three defined transmission rates to a prescribed reference oneβ by using the relative values of the defined one with respect to such a reference transmission rate. In particular, it is easy in that way to characterize the existence of the endemic equilibrium point referring to such a reference transmission rate exceeding a threshold which depends on the remaining parameters including the three relative transmission rates. It will be seen later on in Section 3.2 that the necessary condition β>βc for the reachability of the endemic equilibrium point (equivalent to its existence from its necessary positivity for that concern) is also equivalent to the basic reproduction number to fulfill R0>1, which implies in turn the instability of the disease-free equilibrium point. Thus, one will conclude that if the disease-free equilibrium point is locally asymptotically stable, then the endemic equilibrium one is not reachable. Simulations corroborate that if R0=1, then the disease-free equilibrium point is locally asymptotically stable. As a result, the disease-free equilibrium point is a global attractor if R0≤1 and equivalently if β≤βc.Remark 3.
Note that the sufficient conditionsλ∈0,λc and β>βc of Theorem 2 for the uniqueness of the endemic equilibrium point are independent of the equilibrium subpopulations amount since it is based on the fact that Iee<Nee.
For the particular case,λ=0 which reflects that people, in average, do not go to the doctor under symptoms since the delayed time to go to the doctor is λ−1=+∞. For such a case, it can be deduced as well that the endemic equilibrium is unique. However, the reasoning has to be specific for the discussion of the particular model (4)–(8) obtained by zeroing the terms affected by a factor λ in the differential equations. Since λ=0 yields to division by zero in some relevant equations in the proof of Lemma 1, the next related result is stated and proved separately from Lemma 1 and Theorem 2. Its proof is given in Appendix D.Theorem 3.
Ifλ=0, then there are a unique disease-free equilibrium point and a unique endemic equilibrium. The first one is xdfe=Λ/μ,0,0,0,0T; similarly for λ≠0, and the endemic one satisfies the subsequent relations:(18)Qee=0;Iee=σEeeμ+γ;Ree=γIeeρ+μ=γσρ+μμ+γEee,(19)See=SeeNeeNee=σ+μμ+γΛμβeμ+γ+p2βa+1−p2βsσ,(20)Eee=Λρ+μμ+γZ−σ+μμ+γZμσ+μρ+μ+γσ+μ+ρ,where(21)Z=βeμ+γ+p2βa+1−p2βsσ.
Furthermore, the endemic equilibrium point is reachable if the reference transmission rate exceeds a minimum threshold according to(22)β>βc0=βcλ=0=σ+μ+μ+γμ+γβer+p2βar+1−p2βsrσ.


## 3.2. Local Asymptotic Stability Analysis

The next technical result, proved in Appendix E, relies on the evaluation of the Jacobian matrix with respect to any equilibrium point. It will be then used to characterize the local asymptotic stability of both equilibrium points.Lemma 2.
The Jacobian matrixJ∗=J∗x∗=∂x˜˙t/∂x˜Tt|x∗ of the linearized trajectory x˜t=xt−x∗ around any equilibrium point x∗=S∗,E∗,Q∗,I∗,R∗T is given by J∗=Jij∗∈R5×5 with entries as follows:(23)J11∗=−μ−βeE∗+p2βa+1−p2βsI∗N∗−S∗N∗2,J12∗=−βeS∗N∗−E∗−p2βa+1−p2βsI∗S∗N∗2,J13∗=βeE∗+p2βa+1−p2βsI∗S∗N∗2,J14∗=−p2βa+1−p2βsS∗N∗−I∗−βeE∗S∗N∗2,J15∗=ρ+βeE∗+p2βa+1−p2βsI∗S∗N∗2,J21∗=βeE∗+p2βa+1−p21−λp1βsI∗N∗−S∗N∗2,J22∗=−σ+μ+βeS∗N∗−E∗−p2βa+1−p21−λp1βsI∗S∗N∗2,J23∗=−βeE∗+p2βa+1−p21−λp1βsI∗S∗N∗2,J24∗=p2βa+1−p21−λp1βsS∗N∗−I∗−βeE∗S∗N∗2,J25∗=−βeE∗+p2βa+1−p21−λp1βsI∗S∗N∗2,J31∗=1−p2βsp1λI∗N∗−S∗N∗2,J32∗=−1−p2βsp1λI∗S∗N∗2,J33∗=−μ+α+δ−1−p2βsp1λI∗S∗N∗2,J34∗=λ1−p2+1−p2βsp1λS∗N∗−I∗N∗2,J35∗=−1−p2βsp1λI∗S∗N∗2,J41∗=J43∗=J45∗=0,J42∗=σ,J44∗=−λ1−p2+μ+γ,J51∗=J52∗=0,J53∗=δ,J54∗=γ,J55∗=−ρ+μ.
The above Jacobian matrix can be decomposed asJ∗=T∗−V with the transmission and transition matrices, respectively, −V (being independent of the transmission rates and of the considered equilibrium point) and T∗=J∗+V, where(24)−V=J∗β=0=−μ000ρ0−σ+μ00000−μ+α+δλ1−p200σ0−λ1−p2+μ+γ000δγ−ρ+μ,by using the relative transmission rates βer=βe/β, βsr=βs/β, and βar=βa/β with respect to the reference transmission rate β.The above result is useful for the local stability analysis around the disease-free and endemic equilibrium points, respectively,xdf and xee, via the respective Jacobian matrices Jdf=J∗xdf=T∗xdf−V and Jee=J∗xdf=T∗xee−V. In the calculation of the entries of the Jacobian matrix around x∗, it has been taken into account that the total population is not constant with respect to time due to the mortality associated to the disease. For instance,(25)∂S/N∂Sx=x∗=∂S/S+E+Q+I+R∂Sx=x∗=N∗×1−S∗×1N∗2=1N∗−S∗N∗2,generates one of the additive terms in J11∗ and so on for the remaining entries of the Jacobain matrix around the considered equilibrium point.The following remark gives some considerations of conceptual interest related to the Jacobain matrix around the equilibrium points and its partition into the sum of the transition and transmission matrices.Remark 4.
Note that the transmission matrix−V is a Metzler stability matrix, i.e., its diagonal entries are negative and its off-diagonal entries are nonnegative so that all its eigenvalues are real negative and its associated fundamental (or transition) matrix e−Vt is nonsingular and has nonnegative entries, for all time. Note also that the transmission matrix T∗ is zero if the transmission rates are zero (i.e., it is zero in the absence of disease). The fact that e−Vt has nonnegative entries for all time implies that any solution trajectory is nonnegative for all time under any given nonnegative initial conditions, in the absence of disease which is a necessary condition for the nonnegativity of the whole model to have a nonnegative solution trajectory under nonnegative initial conditions (Theorem 1).
On the other hand,J∗=−VI−V−1T∗=−I−T∗V−1V is stable since −V is stable if and only if the spectral radius of V−1T∗, identical to that of T∗V−1 is small enough. According to equation (23), this fact happens for the disease-free equilibrium point if β is small enough. It can be proved that this is equivalent also to β<βc and it will be seen later on that this is also further equivalent to the basic reproduction number to be less than one. Finally, according to (10), this also implies that the endemic equilibrium point is not reachable.
Local asymptotic stability and instability results around the equilibrium points are stated in the subsequent result, whose proof is given in Appendix F.Theorem 4.
The following properties hold:(i)
The Jacobian matrix around the disease-free equilibrium point is stable if and only if the basic reproduction number(26)R0=ββc−1=βλ1−p2+μ+γβer+p2βar+1−p21−p1λβsrσσ+μλ1−p2+μ+γ<1,with the average critical reference transmission rateβc being defined in (10) and then the disease-free equilibrium point is locally asymptotically stable. Moreover, the spectral radius of V−1Tdf equalizes the basic reproduction number where, such that for the disease-free equilibrium point, the transmission and transition matrices are as follows:(27)Tdf=Tdfβ=Jdf+V=β0−βer0−p2βar+1−p2βsr00βer0λ1−p2p1βsr0000λ1−p21+βsrp100000000000,(28)−V=J∗β=0=−μ000ρ0−σ+μ00000−μ+α+δλ1−p200σ0−λ1−p2+μ+γ000δγ−ρ+μ.(ii)
The disease- free equilibrium point is unstable ifβ>βc⇔R0>1.(iii)
Assume thatλ∈0,λc. Then, the endemic equilibrium point is locally asymptotically stable if and only if β>βc⇔R0>1.(iv)
The above conditions in Properties (ii) and (iii) can be weakened in the sense that the disease-free equilibrium point is locally asymptotically stable if and only ifR0β≤1β≤βc and the endemic one is locally asymptotically stable if and only if R0β≥1β≥βc.Remark 5.
It is now discussed how the calculation of the basic reproduction number can be simplified by manipulating the infective variables, that is, the exposed, infectious, and quarantine subpopulations, only. Note from equations (27)–(28) that the Jacobian matrix at the disease-free equilibrium point is as follows:(29)Jdf=−μ−βe0−p2βa+1−p2βsρ0−σ+μ+βe0p2βa+1−p21−p1λβs000−μ+α+δλ1−p21+βsp100σ0−λ1−p2+μ+γ000δγ−ρ+μ,which makes obvious by its direct inspection that the linearized subsystem of the infective components E, I, and Q around the disease-free equilibrium points is not coupled to the linearized subsystem of noninfective components S and R. Note that the Jacobian matrix of the infective linearized subsystem around the disease-free equilibrium point is given by(30)Jdfi=−VdfiI−Vdfi−1Tdfi=βe−μ+σ0p2βa+1−p21−p1λβs0−μ+α+δλ1−p21+βsp1σ0−λ1−p2+μ+γ,where the corresponding transition and transition matrices are as follows:(31)−Vdfi=−μ+σ000−μ+α+δλ1−p2σ0−λ1−p2+μ+γ;Tdfi=ββer0p2βar+1−p21−p1λβsr00λ1−p2βsrp1000,which concludes that the local stability around the disease-free equilibrium point holds if R0=ρTdfiVdfi<1 which guarantees that Jdfi is a stability matrix since −Vdfi is a stability matrix of eigenvalues −μ+σ, −μ+α+δ, and −λ1−p2+μ+γ. Therefore, the local stability around the disease-free equilibrium point can be analyzed more easily from the Jacobian matrix of the infective linearized subsystem since if Et,It,Qt⟶0 as t⟶∞ for any initial conditions sufficiently close to the disease-free equilibrium point, then the disease-free equilibrium point is locally asymptotically stable and conversely this happens if the eigenvalues of Jdfi have negative real parts.


## 3.3. Global Asymptotic Stability Analysis

The subsequent result extends Theorem4 by proving that the given asymptotic stability results are also global and that only one of the two equilibrium points is a global asymptotically stable attractor depending on the range of values of the basic reproduction number, equivalently, on that of the reference transmission rate for any given fixed triple of relative transmission rates. Its proof is given in Appendix G.Theorem 5.
Assume thatλ∈0,λc. Then, there is only a globally asymptotically stable attractor which is the disease-free equilibrium point if R0β≤1β≤βc and the endemic one if R0>1β>βc.


## 4. Computer Simulations on the Epidemic Evolution

The subsequent parameterization data to test the epidemic model for the COVID-19 pandemic correspond to Spain. It is assumed that the contact ratesce, ca, and cs which conform the respective transmission rates βe, βa, and βs are supposed identical and equal to cr. The transmission probability of the exposed pe is 10−3 times lower than that of the symptomatic ps which equalizes that of the asymptomatic individuals pa. The average time before an infectious individual goes to the doctor is λ−1=2 days. The tracing effectiveness p1 is 0.5 while the average time of staying in quarantine is equal to δ−1=σ−1+γ−1. The remaining model parameters and the initial conditions are displayed in Tables 1 and 2.Table 1
Values of the model parameters.
ParameterDefinitionValueCiteΛ/NBirth rate2.088× 10−5 day−1[9]crContact rate10.58 day−1[23]psTransmission probability of infected people0.23[29]p2Probability to be asymptomatic0.16[2, 13]μNatural death rate2.282× 10−5 day−1[9]1/σIncubation period7 days[23]γRemoval rate0.1 day−1[3]αDisease-induced death rate0.01 day−1[23]Table 2
Initial conditions.
VariableDefinitionValueS0Susceptible population46× 106E0Exposed population100Q0Quarantine population1I0Infected population1R0Recovered population0The obtained simulation results are depicted in Figure2(a); it shows the dynamics of the variables S, R, Q, E, and I and the cumulative number of deaths with respect to time. The variables show characteristics of an underdamped dynamical system; its response is a sum of a transient response (which is damped sinusoidal) and the steady state. In this case, the disease persists, which is an expected result since R0 = 2.36. It is important to note that, approximately, 300 days after the outbreak, the total number of deaths reaches approximately 1 million, and at the end of the simulation (600 days after the outbreak) is doubled. This high mortality can be attributed to the high number of contagious contacts. The disease spreading can be decreased by reducing the value of cr=ce=cs=ca; therefore, in an alternative simulation, cr was reduced to 4 (R0 = 0.89), and as it can be seen in Figure 2(b), the disease tends to disappear; the variables Q, E, and I show an overdamped response, and their values continuously decrease with respect to time. As long as their values and the total number of deaths are small compared with S long over the simulation, the variable S seems to be unchanging (see the top graph of Figure 2(b)). The results of both above simulations correspond to the expected behavior in the cases of which the basic reproduction number is greater and lower than unity. Besides, considering the cumulative number of deaths of both simulations, it is clear that reducing the value of the parameter of average contact rate cr could prevent a substantial increase in the total deaths.Figure 2
S, R, E, I, and Q population sizes with respect to time when cr = 10.6 and R0 = 2.360 (a) and cr = 4 and R0 = 0.89 (b).
(a)(b)On the other hand, Figure3 shows the disease evolution as the basic reproduction number is unity of close to unity. It is seen that if the reproduction number is unity so that the disease-free and the endemic equilibrium points coincide, the disease-free is asymptotically stable as in the case when it takes values below unity as expected from the theoretical discussion performed in the former section. It is also seen that, for values of the basic reproduction number which are slightly higher than unity, the endemic equilibrium is a globally asymptotically stable attractor.Figure 3
Comparison of the infection evolution for the basic reproduction number close to unity: (a) infection evolution for the unity basic reproduction number, (b) infection evolution for the basic reproduction number slightly below unity, and (c) infection evolution for the basic reproduction number slightly over unity.
(a)(b)(c)


## 5. A Case Study: Evaluation of the Model with Data on COVID-19 Taken from Spanish Cantabria Community

### 5.1. Data Acquisition

This section evaluates the proposed model with real discrete data from Cantabria Community. The tuple of daily discrete data of thek−th day is defined yk=tk,pk,dk,rk where tk is the number of performed tests, pk is the number of positive tested results, dk is the number of deaths, and rk is the number of recovered individuals. Over the above array, the subindices d and c are used for daily and cumulative data. From cumulative data, one can obtain the daily ones by using the derivative of the accumulated ones by Euler forward approximation which results in ydk+h=yck+h−ych/h, and typically h will be taken as a period of one day. Some relevant data from Cantabria, from February 29, 2020, to February 14, 2021, are displayed in Figure 4.Figure 4
Epidemiological situation in Cantabria.It can be pointed out that it is of interest to evaluate from the data the number of individuals in quarantine, which also includes the hospitalized ones since, approximately, the deaths are proportional to the quarantined individuals, i.e.,dk=αQk and the fraction of recovered individuals is proportional to the quarantined individuals, i.e., rk=δQk. Thus, by discretizing (8) for one day period, one gets(32)rk=δQk=Rk+1−Rk+ρ+μRk−γIk.In the first propagation wave of the disease, since the number of infectiousIk individuals compared with recovered individuals Rk is close to zero and ρ+μ≪1, it can be stated that rk is approximately the daily increment of the recovered since from the above relation, rk=δQk≅Rk+1−Rk. Note that rk is a part of the recovered individuals, that is, those recovered individuals from quarantine, since most of the recovered individuals from the slight or the asymptomatic infection are not accounted for the data processing. Therefore, we keep in the following lowercase notation rk for this partial account of the estimated recovered daily increment Rk+1−Rk. The period from February 29, 2020, to February 14, 2021, gave a cumulative of recovered and death people which was 25,500 which is only a 5% of the population of Cantabria. Therefore, S/N has been taken unity along that period. Regarding the number of positive tests, death individuals, and recovered tested daily individuals pk, dk, and rk, some average delays τ1,τ2, and τ3 exist from the contributions of the daily infectious and quarantined individuals caused by the transmission process and the delays in data processing. The used evolution equations are as follows:(33)pk=p^k+wpk=λ1−p21+p1pscrIk−τ1+wpk,(34)dk=d^k+wdk=α+μQk−τ2+wdk,(35)rk=r^k+wrk=δQk−τ3+wrk,where the superscript hat stands for the various estimates and wpk0∞, wdk0∞, and wrk0∞ are disturbance white noises. To reduce noise effects, the centered moving-average filter which takes 14 samples has been applied to (33)–(35). By combining the above equations, one gets for the discretization of the estimate of the quarantined individuals through (6) for h=1 via the Euler forward approximation that(36)Q^k+1=Q^k+λ1−p21+p1pscrIk−α+μQ^k−δQ^k=Q^k+λ1−p21+p1pscrpk+τ1−wpk+τ1λ1−p21+p1pscr−α+μα+μdk+τ2−wdk+τ2+wrk+τ3−rk+τ3=Q^k+pk+τ1−dk+τ2−rk+τ3+wdk+τ2+wrk+τ3−wpk+τ1=Q^k+p^k+τ1−d^k+τ2−r^k+τ3.


### 5.2. Estimation of the Delays

Note that as far as the coexistence with COVID-19 increases, it is expected faster responses in the performance of positivity tests and public supply of data. Therefore, the potential estimation of delays can be adapted to various disease evolution periods. Now, defineτ4=τ2−τ3 and τ5=τ1−τ2. Algorithm 1 has been performed to estimate them from the registered data.Algorithm 1.
(estimation ofτ4=τ2−τ3 and τ5=τ1−τ2)Step 1. Initialization: Fix a natural number n and a set T5 of “a priori” potential delays τ5 as T5=−n,−n+1,⋯,n−1,n and make ℓ←1 and τℓ5i=−n for i=1,2.Divide the data setpk,dk,rk into subsets for different disease evolution time periods. In particular, p1k,d1k,r1k for the period from February 29, 2020, to July 17, 2020, and to estimate τ41,τ51, and p2k,d2k,r2k for the period from July 18, 2020, to December 31, 2021, to estimate τ42,τ52.Step 2. For each subset, seek for delays between r^. and d^. (in particular, the delays between the corresponding additive contributions of d^k+τ2 and r^k+τ3 to (45)) with MATLAB finddelay() (for delays between signals) to get τ4i for i=1,2.Step 3. Compute τ6i, i=1,2, with finddelay() to align Q^k−τ2i with respect to d^k and r^k−τ4i since they are both proportional to Q^k−τ2i. Therefore, τ6i=finddelayQ^k−τ2i,d^k=finddelayQ^k−τ2i,r^k−τ4i, i=1,2.Step 4. Calculate p^k−τ5i from (33).Step 5. Calculate (36) under the equivalent form as follows:(37)Q^k+1−τ2i=Q^k−τ2i+p^k+τ5i−d^k−r^k−τ4i.Step 6. Apply the linear regression method to the data sets Q^k−τ2i,r^k−τ4i and Q^k−τ2i,d^k and compute the coefficients of determination Rℓi2r and Rℓi2d, i=1,2.Step 7. Compute fℓi=Rℓi2r+Rℓi2d/2.Step 8. While ℓ≤2n do ℓ←ℓ+1 and go to Step 4.Step 9. τ5i=Argmaxfℓiτℓ5i∈T5.Step 10. End.Remark 6.
Both coefficients of determinationRℓi2r and Rℓi2d, i=1,2, are the proportions of variability in the set of given data due to prediction [35, 36], and they range from zero to unity. The data become better fitted as the corresponding coefficient of determination becomes close to unity. In contrast, values close to zero lead to a high variability and a bad prediction as a result.
As a result of applying Algorithm1, one obtains τj=τji;j=4,5,6,i=1,2, with the values τ41=26days, τ51=11days, and τ61=9days and τ42=0days, τ52=−1days, and τ62=13days, and the number of individuals in quarantine due to positive tests is displayed in Figure 5.Figure 5
Estimated people in quarantine after testing positive for COVID-19 at different time periods.


### 5.3. Estimation of Average Mortality Rate and Stay in Quarantine Periods

The estimations ofαi and δi, i=1,2, are displayed in Table 3. Note that α1>α2 as expected from a higher hospitalization percentage along the first studied period. Regarding the recovery rate, after the first wave, tracing techniques were implemented when someone with symptoms went to hospital and his/her close contacts were tested to detect new infections, and those tested positive were quarantined which can explain that δ2<δ1.Table 3
Estimated parameters via linear regression.
ParameterValueErrorR2α10.00630.00010.96α27.1 ∗ 10−40.3 ∗ 10−40.75δ10.0570.8 ∗ 10−40.97δ20.0430.0020.80Figures6 and 7 display the linear regressions along the first and second evaluation time periods to estimate the mortality rates and stay in quarantine period.Figure 6
Linear regression in the period from 01/03/2020 to 16/07/2020 to estimate the parametersα1 and δ1 values.Figure 7
Linear regression in the period from 17/07/2020 to 02/02/2021 to estimate the parametersα2 and δ2.


### 5.4. Estimations ofλ,cr, and p1

By approximatingS/N to unity, we take equations (5) and (7), the infective subsystem of the epidemic model, and we rewrite the set of third-order equations of first-order as a differential equation of third-order with the quarantined subpopulation as variable. This leads to(38)Q…t+b1Q¨t+b2Q˙t−b3Qt=0,where(39)b1=λ1−p2+γ+α+δ+σ+3μ−crpe,b2=α+δ+μλ1−p2+γ+μ+σ+μ−crpeλ1−p2+γ+α+δ+2μ−σcrpsp2+1−p21−p1λ,b3=crpe−σ−μλ1−p2+γ+μ+σcrpsp2+1−p21−p1λα+δ+μ,where pe=βe/cr and it has been considered that pa=ps=βs/cr depending on the average transmission rates βe=βa and βs and that the average contact rates satisfy cr=ce=cs=ca. It has been also assumed that λ is the average rate of visits to the doctor or to hospital for the particular disease COVID-19. Now, apply the Euler forward approximation to (38) for discretization which leads to(40)Qk=3−b1Qk−1+2b1−b2−3Qk−2+1−b1+b2+b3Qk−3,and equalize the above expression by equating corresponding coefficients to the prediction values of the quarantined individuals given by(41)Q^k=−a^1Qk−1−a^2Qk−1−a^Q3k−3,which leads to the system of algebraic equations J1x^1=v^1 with x^1=λ^,c^r,λ^c^rT which are in expanded form:(42)1−p2−pe0−σ+μ1−p2peγ+μ+psσ1−p2pe−psp1σα+δ+σ+2μ1−p2−peγ+α+δ+2μ+psσ1−p2psp1σ−peλ^c^rλ^c^r=a^1+3−γ+ε+δ+σ+3μ−a^1+a^2+a^3α+δ+μ+σ+μγ+μ2a^1+a^2+3−α+δγ−σγ+α+δ−μ2γ+2α+2δ+2σ+3μ,where the parameters which are not unknowns are taken from Table 1. Thus, since the third row is not independent of the two former ones, one takes the minimum square error solution λ^,c^r of (42) given by(43)λ^,c^r=Argminp1∈R+J1x^1−v^12,with stopping computational value error less than 10−4. Now, proceed in the same way to get a second-order differential equation of It, by using the combination of (5) and (7), resulting in(44)I¨t=c1I˙t+c2It,with(45)c1=crp−eσ−λ1−p2−γ−2μ,(46)c2=σcrp2pe+1−p21−λp1ps+crpe−σ+μλ1−p2+μ+γ.By applying the Euler forward discretization method, one gets(47)Ik=2+c1Ik−1+c2−c1−1Ik−2.In this case, one gets proceeding in the same way as done with (40) to get (42) the algebraic system J2p^1=v^2 with the equivalent explicit expanded scalar relationships as follows:(48)2+c1=−λ^1−p21+c^rpsp1d^1,(49)c2−c1−1=−λ^1−p21+c^rpsp1d^2.And the estimationp^1 of the effectiveness tracing p1 is calculated by replacing the necessary values of Table 1 and the solution λ^,c^r of (43) results to be(50)p^1=Argminp1∈R+J2x^1−v^22,with stopping computational error less than 10−4. The solution is found through the Newton–Raphson method [42, 43].


### 5.5. Relationship between Positive and Total Tests Interpreted via the Estimations ofcr and p1

The total and positive tests are described, respectively, by(51)tk=1+crkp1kSkNkλ1−p2Ik+εk,(52)pk=1+pscrkp1kSkNkλ1−p2Ik.As a consequence of testing symptomatic people with COVID-191−p2λIk, we traced contacts crλp1k1−p2IS/N and other people that present similar symptoms to COVID-19 but are not infected εk. By taking Sk/Nk≅1, one gets(53)tk=rakpk+εk,where(54)rak=1+crkp1k1+crkp1kps.The estimations of the values ofλ, cr, and p1 for each one of the mentioned time intervals are displayed in Table 4. The smaller values of the estimation of λ were found in the period from 19/12/2020 to 06/01/2021 which matches with winter holidays when a significant number of people moved to their relative houses so that people under infection symptoms have found more difficulty to go to the outpatient. Regarding the average number of contacts per day cr, as tighter restrictions were applied, smaller values of cr were observed. The smallest value c^r3=2.07 is found from 19/03/2020 to 28/05/2020 approximately as the quarantine was imposed while the next closer smaller one c^r12=3.14 was reached from 25/11/2020 to 13/12/2020 as the state of alertness was still valid. Thus, these values are closely related to the levels of imposed restrictions.Table 4
Estimated parametrical values.
DatesNatural positivity rateλ (days−1)Average number of close contacts per daycrTracing the effectivenessp1[29/02/2020, 10/03/2020]λ^1 = 0.685c^r1 = 12.88p^11 = 0.06[11/03/2020, 18/03/2020]λ^2 = 0.643c^r2 = 4.88p^12 = 0.9[19/03/2020, 28/05/2020]λ^3 = 0.938c^r3 = 2.07p^13 = 0.25[28/05/2020, 30/06/2020]λ^4 = 0.859c^r4 = 4.15p^14 = 0.9[17/07/2020, 06/08/2020]λ^5 = 0.930c^r5 = 7.55p^15 = 0.06[10/08/2020, 02/09/2020]λ^6 = 0.932c^r6 = 4.81p^16 = 0.02[04/09/2020, 17/09/2020]λ^7 = 0.737c^r7 = 4.63p^17 = 0.42[17/09/2020, 14/10/2020]λ^8 = 0.871c^r8 = 3.21p^18 = 0.28[15/10/2020, 30/10/2020]λ^9 = 0.932c^r9 = 6.10p^19 = 0.14[30/10/2020, 09/11/2020]λ^10 = 0.933c^r10 = 5.57p^110 = 0.25[09/11/2020, 23/11/2020]λ^11 = 0.888c^r11 = 4.68p^111 = 0.36[25/11/2020, 13/12/2020]λ^12 = 0.744c^r12 = 3.14p^112 = 1[19/12/2020, 06/01/2021]λ^13 = 0.628c^r13 = 6.38p^113 = 0.78[07/01/2021, 20/01/2021]λ^14 = 0.881c^r14 = 5.06p^114 = 0.05Note that the parameterλ−1 of time of meeting doctor or health service for testing varies approximately from 1.58 days to 1.07 days. Assuming that, at first, COVID-19 was underestimated and its symptoms were mixed up with cold/flu ones, infected people with symptoms might typically have to visit the doctors at lower rates. On the other hand, holiday seasons could also be associated with lower values of λ since people might be outside their usual homes so that health assistance might be reduced. The tracing effectiveness p1 was never before implemented in Spain on big scale cases such as in this pandemic; therefore, there has been notably uncertainty regarding how to implement it. In that way, as the tested positivity was high, people were tested randomly. The variances of the parameters λ^k, c^rk, and p^1k are calculated through the values of Table 4. On the other hand, Figures 8–10 display the quarantined, death, and recovered individuals in both periods compared with their estimations taking account of the variances of the involved parameters, for instance, and concerning Figure 8, the corresponding variances, or its associated typical deviations, of Q^k are obtained via (36) as follows. For each k−th sample, corresponding to a particular day, one takes data from two weeks (that is, 14 days) neighboring days, namely, the former seven ones, the current ones, and the next six ones so there is one delay of six days to provide the evaluation of the mean value of the involved data over each day. One proceeds in this way for the estimations on the mentioned groups of two weeks for the estimations of dk, rk, and pk. With those data, one calculates the k-th variance, or its associate typical deviation, of the estimation of the quarantined subpopulation Qk.Figure 8
Comparison between measured people in quarantine per day (black dots) with respect to the estimated people in quarantine per day (blue solid line) with the estimated parameters. The dashed red lines show the standard error of the simulation.Figure 9
Comparison between measured death people per day (black dots) with respect to the estimation (blue solid line) carried out with the estimated parameters. The dashed red lines show the standard error of the simulation.Figure 10
Comparison between recovered people per day (black dots) with respect to the estimation (blue solid line) carried out with the estimated parameters. The dashed red lines show the standard error of the solution.Figure11 shows a linear regression of total tests tk versus positive tests pk in two periods of time from 01/03/2020 to16/07/2020 and from 17/07/2020 to 02/02/2021. The adjustment of the corresponding linear regression equation tk=rakpk+εk via least-squares gives the pairs rak,εk. Note from (54) that rak∈0,1 and that tk=rakpk+εk≥pk so that εk≥1−rakpak≥0 so that, for potentially got negative values of rak and εk, the corresponding data are not collected.Figure 11
Linear regression of total tests versus positive tests in the periods of time from 01/03/2020 to 16/07/2020 (a) and from 17/07/2020 to 02/02/2021 (b).
(a)(b)


## 5.1. Data Acquisition

This section evaluates the proposed model with real discrete data from Cantabria Community. The tuple of daily discrete data of thek−th day is defined yk=tk,pk,dk,rk where tk is the number of performed tests, pk is the number of positive tested results, dk is the number of deaths, and rk is the number of recovered individuals. Over the above array, the subindices d and c are used for daily and cumulative data. From cumulative data, one can obtain the daily ones by using the derivative of the accumulated ones by Euler forward approximation which results in ydk+h=yck+h−ych/h, and typically h will be taken as a period of one day. Some relevant data from Cantabria, from February 29, 2020, to February 14, 2021, are displayed in Figure 4.Figure 4
Epidemiological situation in Cantabria.It can be pointed out that it is of interest to evaluate from the data the number of individuals in quarantine, which also includes the hospitalized ones since, approximately, the deaths are proportional to the quarantined individuals, i.e.,dk=αQk and the fraction of recovered individuals is proportional to the quarantined individuals, i.e., rk=δQk. Thus, by discretizing (8) for one day period, one gets(32)rk=δQk=Rk+1−Rk+ρ+μRk−γIk.In the first propagation wave of the disease, since the number of infectiousIk individuals compared with recovered individuals Rk is close to zero and ρ+μ≪1, it can be stated that rk is approximately the daily increment of the recovered since from the above relation, rk=δQk≅Rk+1−Rk. Note that rk is a part of the recovered individuals, that is, those recovered individuals from quarantine, since most of the recovered individuals from the slight or the asymptomatic infection are not accounted for the data processing. Therefore, we keep in the following lowercase notation rk for this partial account of the estimated recovered daily increment Rk+1−Rk. The period from February 29, 2020, to February 14, 2021, gave a cumulative of recovered and death people which was 25,500 which is only a 5% of the population of Cantabria. Therefore, S/N has been taken unity along that period. Regarding the number of positive tests, death individuals, and recovered tested daily individuals pk, dk, and rk, some average delays τ1,τ2, and τ3 exist from the contributions of the daily infectious and quarantined individuals caused by the transmission process and the delays in data processing. The used evolution equations are as follows:(33)pk=p^k+wpk=λ1−p21+p1pscrIk−τ1+wpk,(34)dk=d^k+wdk=α+μQk−τ2+wdk,(35)rk=r^k+wrk=δQk−τ3+wrk,where the superscript hat stands for the various estimates and wpk0∞, wdk0∞, and wrk0∞ are disturbance white noises. To reduce noise effects, the centered moving-average filter which takes 14 samples has been applied to (33)–(35). By combining the above equations, one gets for the discretization of the estimate of the quarantined individuals through (6) for h=1 via the Euler forward approximation that(36)Q^k+1=Q^k+λ1−p21+p1pscrIk−α+μQ^k−δQ^k=Q^k+λ1−p21+p1pscrpk+τ1−wpk+τ1λ1−p21+p1pscr−α+μα+μdk+τ2−wdk+τ2+wrk+τ3−rk+τ3=Q^k+pk+τ1−dk+τ2−rk+τ3+wdk+τ2+wrk+τ3−wpk+τ1=Q^k+p^k+τ1−d^k+τ2−r^k+τ3.


## 5.2. Estimation of the Delays

Note that as far as the coexistence with COVID-19 increases, it is expected faster responses in the performance of positivity tests and public supply of data. Therefore, the potential estimation of delays can be adapted to various disease evolution periods. Now, defineτ4=τ2−τ3 and τ5=τ1−τ2. Algorithm 1 has been performed to estimate them from the registered data.Algorithm 1.
(estimation ofτ4=τ2−τ3 and τ5=τ1−τ2)Step 1. Initialization: Fix a natural number n and a set T5 of “a priori” potential delays τ5 as T5=−n,−n+1,⋯,n−1,n and make ℓ←1 and τℓ5i=−n for i=1,2.Divide the data setpk,dk,rk into subsets for different disease evolution time periods. In particular, p1k,d1k,r1k for the period from February 29, 2020, to July 17, 2020, and to estimate τ41,τ51, and p2k,d2k,r2k for the period from July 18, 2020, to December 31, 2021, to estimate τ42,τ52.Step 2. For each subset, seek for delays between r^. and d^. (in particular, the delays between the corresponding additive contributions of d^k+τ2 and r^k+τ3 to (45)) with MATLAB finddelay() (for delays between signals) to get τ4i for i=1,2.Step 3. Compute τ6i, i=1,2, with finddelay() to align Q^k−τ2i with respect to d^k and r^k−τ4i since they are both proportional to Q^k−τ2i. Therefore, τ6i=finddelayQ^k−τ2i,d^k=finddelayQ^k−τ2i,r^k−τ4i, i=1,2.Step 4. Calculate p^k−τ5i from (33).Step 5. Calculate (36) under the equivalent form as follows:(37)Q^k+1−τ2i=Q^k−τ2i+p^k+τ5i−d^k−r^k−τ4i.Step 6. Apply the linear regression method to the data sets Q^k−τ2i,r^k−τ4i and Q^k−τ2i,d^k and compute the coefficients of determination Rℓi2r and Rℓi2d, i=1,2.Step 7. Compute fℓi=Rℓi2r+Rℓi2d/2.Step 8. While ℓ≤2n do ℓ←ℓ+1 and go to Step 4.Step 9. τ5i=Argmaxfℓiτℓ5i∈T5.Step 10. End.Remark 6.
Both coefficients of determinationRℓi2r and Rℓi2d, i=1,2, are the proportions of variability in the set of given data due to prediction [35, 36], and they range from zero to unity. The data become better fitted as the corresponding coefficient of determination becomes close to unity. In contrast, values close to zero lead to a high variability and a bad prediction as a result.
As a result of applying Algorithm1, one obtains τj=τji;j=4,5,6,i=1,2, with the values τ41=26days, τ51=11days, and τ61=9days and τ42=0days, τ52=−1days, and τ62=13days, and the number of individuals in quarantine due to positive tests is displayed in Figure 5.Figure 5
Estimated people in quarantine after testing positive for COVID-19 at different time periods.


## 5.3. Estimation of Average Mortality Rate and Stay in Quarantine Periods

The estimations ofαi and δi, i=1,2, are displayed in Table 3. Note that α1>α2 as expected from a higher hospitalization percentage along the first studied period. Regarding the recovery rate, after the first wave, tracing techniques were implemented when someone with symptoms went to hospital and his/her close contacts were tested to detect new infections, and those tested positive were quarantined which can explain that δ2<δ1.Table 3
Estimated parameters via linear regression.
ParameterValueErrorR2α10.00630.00010.96α27.1 ∗ 10−40.3 ∗ 10−40.75δ10.0570.8 ∗ 10−40.97δ20.0430.0020.80Figures6 and 7 display the linear regressions along the first and second evaluation time periods to estimate the mortality rates and stay in quarantine period.Figure 6
Linear regression in the period from 01/03/2020 to 16/07/2020 to estimate the parametersα1 and δ1 values.Figure 7
Linear regression in the period from 17/07/2020 to 02/02/2021 to estimate the parametersα2 and δ2.


## 5.4. Estimations ofλ,cr, and p1

By approximatingS/N to unity, we take equations (5) and (7), the infective subsystem of the epidemic model, and we rewrite the set of third-order equations of first-order as a differential equation of third-order with the quarantined subpopulation as variable. This leads to(38)Q…t+b1Q¨t+b2Q˙t−b3Qt=0,where(39)b1=λ1−p2+γ+α+δ+σ+3μ−crpe,b2=α+δ+μλ1−p2+γ+μ+σ+μ−crpeλ1−p2+γ+α+δ+2μ−σcrpsp2+1−p21−p1λ,b3=crpe−σ−μλ1−p2+γ+μ+σcrpsp2+1−p21−p1λα+δ+μ,where pe=βe/cr and it has been considered that pa=ps=βs/cr depending on the average transmission rates βe=βa and βs and that the average contact rates satisfy cr=ce=cs=ca. It has been also assumed that λ is the average rate of visits to the doctor or to hospital for the particular disease COVID-19. Now, apply the Euler forward approximation to (38) for discretization which leads to(40)Qk=3−b1Qk−1+2b1−b2−3Qk−2+1−b1+b2+b3Qk−3,and equalize the above expression by equating corresponding coefficients to the prediction values of the quarantined individuals given by(41)Q^k=−a^1Qk−1−a^2Qk−1−a^Q3k−3,which leads to the system of algebraic equations J1x^1=v^1 with x^1=λ^,c^r,λ^c^rT which are in expanded form:(42)1−p2−pe0−σ+μ1−p2peγ+μ+psσ1−p2pe−psp1σα+δ+σ+2μ1−p2−peγ+α+δ+2μ+psσ1−p2psp1σ−peλ^c^rλ^c^r=a^1+3−γ+ε+δ+σ+3μ−a^1+a^2+a^3α+δ+μ+σ+μγ+μ2a^1+a^2+3−α+δγ−σγ+α+δ−μ2γ+2α+2δ+2σ+3μ,where the parameters which are not unknowns are taken from Table 1. Thus, since the third row is not independent of the two former ones, one takes the minimum square error solution λ^,c^r of (42) given by(43)λ^,c^r=Argminp1∈R+J1x^1−v^12,with stopping computational value error less than 10−4. Now, proceed in the same way to get a second-order differential equation of It, by using the combination of (5) and (7), resulting in(44)I¨t=c1I˙t+c2It,with(45)c1=crp−eσ−λ1−p2−γ−2μ,(46)c2=σcrp2pe+1−p21−λp1ps+crpe−σ+μλ1−p2+μ+γ.By applying the Euler forward discretization method, one gets(47)Ik=2+c1Ik−1+c2−c1−1Ik−2.In this case, one gets proceeding in the same way as done with (40) to get (42) the algebraic system J2p^1=v^2 with the equivalent explicit expanded scalar relationships as follows:(48)2+c1=−λ^1−p21+c^rpsp1d^1,(49)c2−c1−1=−λ^1−p21+c^rpsp1d^2.And the estimationp^1 of the effectiveness tracing p1 is calculated by replacing the necessary values of Table 1 and the solution λ^,c^r of (43) results to be(50)p^1=Argminp1∈R+J2x^1−v^22,with stopping computational error less than 10−4. The solution is found through the Newton–Raphson method [42, 43].


## 5.5. Relationship between Positive and Total Tests Interpreted via the Estimations ofcr and p1

The total and positive tests are described, respectively, by(51)tk=1+crkp1kSkNkλ1−p2Ik+εk,(52)pk=1+pscrkp1kSkNkλ1−p2Ik.As a consequence of testing symptomatic people with COVID-191−p2λIk, we traced contacts crλp1k1−p2IS/N and other people that present similar symptoms to COVID-19 but are not infected εk. By taking Sk/Nk≅1, one gets(53)tk=rakpk+εk,where(54)rak=1+crkp1k1+crkp1kps.The estimations of the values ofλ, cr, and p1 for each one of the mentioned time intervals are displayed in Table 4. The smaller values of the estimation of λ were found in the period from 19/12/2020 to 06/01/2021 which matches with winter holidays when a significant number of people moved to their relative houses so that people under infection symptoms have found more difficulty to go to the outpatient. Regarding the average number of contacts per day cr, as tighter restrictions were applied, smaller values of cr were observed. The smallest value c^r3=2.07 is found from 19/03/2020 to 28/05/2020 approximately as the quarantine was imposed while the next closer smaller one c^r12=3.14 was reached from 25/11/2020 to 13/12/2020 as the state of alertness was still valid. Thus, these values are closely related to the levels of imposed restrictions.Table 4
Estimated parametrical values.
DatesNatural positivity rateλ (days−1)Average number of close contacts per daycrTracing the effectivenessp1[29/02/2020, 10/03/2020]λ^1 = 0.685c^r1 = 12.88p^11 = 0.06[11/03/2020, 18/03/2020]λ^2 = 0.643c^r2 = 4.88p^12 = 0.9[19/03/2020, 28/05/2020]λ^3 = 0.938c^r3 = 2.07p^13 = 0.25[28/05/2020, 30/06/2020]λ^4 = 0.859c^r4 = 4.15p^14 = 0.9[17/07/2020, 06/08/2020]λ^5 = 0.930c^r5 = 7.55p^15 = 0.06[10/08/2020, 02/09/2020]λ^6 = 0.932c^r6 = 4.81p^16 = 0.02[04/09/2020, 17/09/2020]λ^7 = 0.737c^r7 = 4.63p^17 = 0.42[17/09/2020, 14/10/2020]λ^8 = 0.871c^r8 = 3.21p^18 = 0.28[15/10/2020, 30/10/2020]λ^9 = 0.932c^r9 = 6.10p^19 = 0.14[30/10/2020, 09/11/2020]λ^10 = 0.933c^r10 = 5.57p^110 = 0.25[09/11/2020, 23/11/2020]λ^11 = 0.888c^r11 = 4.68p^111 = 0.36[25/11/2020, 13/12/2020]λ^12 = 0.744c^r12 = 3.14p^112 = 1[19/12/2020, 06/01/2021]λ^13 = 0.628c^r13 = 6.38p^113 = 0.78[07/01/2021, 20/01/2021]λ^14 = 0.881c^r14 = 5.06p^114 = 0.05Note that the parameterλ−1 of time of meeting doctor or health service for testing varies approximately from 1.58 days to 1.07 days. Assuming that, at first, COVID-19 was underestimated and its symptoms were mixed up with cold/flu ones, infected people with symptoms might typically have to visit the doctors at lower rates. On the other hand, holiday seasons could also be associated with lower values of λ since people might be outside their usual homes so that health assistance might be reduced. The tracing effectiveness p1 was never before implemented in Spain on big scale cases such as in this pandemic; therefore, there has been notably uncertainty regarding how to implement it. In that way, as the tested positivity was high, people were tested randomly. The variances of the parameters λ^k, c^rk, and p^1k are calculated through the values of Table 4. On the other hand, Figures 8–10 display the quarantined, death, and recovered individuals in both periods compared with their estimations taking account of the variances of the involved parameters, for instance, and concerning Figure 8, the corresponding variances, or its associated typical deviations, of Q^k are obtained via (36) as follows. For each k−th sample, corresponding to a particular day, one takes data from two weeks (that is, 14 days) neighboring days, namely, the former seven ones, the current ones, and the next six ones so there is one delay of six days to provide the evaluation of the mean value of the involved data over each day. One proceeds in this way for the estimations on the mentioned groups of two weeks for the estimations of dk, rk, and pk. With those data, one calculates the k-th variance, or its associate typical deviation, of the estimation of the quarantined subpopulation Qk.Figure 8
Comparison between measured people in quarantine per day (black dots) with respect to the estimated people in quarantine per day (blue solid line) with the estimated parameters. The dashed red lines show the standard error of the simulation.Figure 9
Comparison between measured death people per day (black dots) with respect to the estimation (blue solid line) carried out with the estimated parameters. The dashed red lines show the standard error of the simulation.Figure 10
Comparison between recovered people per day (black dots) with respect to the estimation (blue solid line) carried out with the estimated parameters. The dashed red lines show the standard error of the solution.Figure11 shows a linear regression of total tests tk versus positive tests pk in two periods of time from 01/03/2020 to16/07/2020 and from 17/07/2020 to 02/02/2021. The adjustment of the corresponding linear regression equation tk=rakpk+εk via least-squares gives the pairs rak,εk. Note from (54) that rak∈0,1 and that tk=rakpk+εk≥pk so that εk≥1−rakpak≥0 so that, for potentially got negative values of rak and εk, the corresponding data are not collected.Figure 11
Linear regression of total tests versus positive tests in the periods of time from 01/03/2020 to 16/07/2020 (a) and from 17/07/2020 to 02/02/2021 (b).
(a)(b)


## 6. Conclusions

This paper has proposed a SEQIR model which incorporates a quarantined subpopulation to those of the standard ones in the common SEIR-type models. The quarantined subpopulation contains the individuals who have been tested positive. This includes the hospitalized individuals including those under intensive care and also those quarantined at their homes having no serious symptoms. The model parameterizes the average transmission depending on the infective contacts as the average rate of attendance to the doctor for disease positivity testing or symptom evaluation checking and the average effectiveness tracing on positive cases.The proposed model considers that the proportions of asymptomatic and symptomatic individuals in the quarantined and the infectious subpopulations can potentially have different average transmission rates of contagion to the susceptible group. Also, various average transmission rates are assumed to be parameterized by the infective contacts between individuals which allows to evaluate the infection transmission reduction under public intervention measures like confinements or mobility restrictions, social distance maintenance, or use of masks. A main novelty of the proposed model is that it is parameterized by the average rate of attendance to the doctor for disease positivity testing or evaluation of symptoms checking and by the effectiveness in tracing the positive cases. It is proved that the increase in the attendance rate for testing of symptomatic suspect individuals contributes positively to the incorporation to the quarantined subpopulation of those being positively tested so that the infection transmission might be controlled more efficiently.The model is also proved to be positive, an obvious need for the coherency of any epidemic model, in the sense that the state trajectory solution is seen to be nonnegative for all time for any given set of nonnegative initial conditions while it is proved to have unique disease-free and endemic equilibrium points. The disease-free equilibrium point is proved to be locally asymptotically stable if the basic reproduction number is less than unity, and under these conditions, the endemic one is not reachable, while it is locally asymptotically stable if such a value exceeds unity. Both points are globally asymptotically stable under the abovementioned conditions, and they are locally asymptotically stable. The mathematical condition for the basic reproduction number to be less than or larger than unity is proved to be equivalent to the transmission rates to be less than or larger than a certain critical transmission rate under the assumption that the remaining model parameters are kept unmodified. This helps to interpret the eventual transmission attenuation in terms of reduction of the disease transmission rates, implying as a result a decrease in the infective contagion contacts, as being dependent on the potential intervention measures. The model is tested by computer simulations for COVID-19 with data from Spain. Subsequently, a detailed case study was performed against COVID-19 through official data taken from the Spanish Cantabria Autonomous Community. In particular, the estimations of the average periods lasted by people to get to a health center or to a doctor for eventual positivity inspection, the positivity efficiency of tests, and the average number of contagion contacts are estimated from real data under the confinement period and also for a subsequent period. In this second period, the first strong public intervention measures of almost total general confinement were left, while still the population was under severe rules like general use of masks, social distance keeping, closing of nocturnal amusement places, or very important limitation rules on the number of attendees to spectacles and restaurants.In the proposed model, there is a good trade-off between accuracy, transparency, and flexibility. Even if the quarantined subpopulation gathers people who were tested positive and were staying either at home or in hospital, reducing so the model complexity, a quantitative and qualitative fit of data, with respect to real data, has been checked. Regarding the transparency, the calculated basic reproduction number shows a simple relation with respect to parameters such as the average number of close contacts per day or the tracing effectiveness. The proposed model has been then evaluated with real data acquisition from the Spanish Autonomous Region of Cantabria. From the data source, the number of positive test result per day and the number of death people due to the disease and recovered people were gathered.A noteworthy delay between the quarantined with respect to death and recovered can be interpreted in terms of delays in providing data and in the eventual attendance for checking eventual positivity testing. The analysis of the imposed restriction with respect to time gives an overview of how the limitations might affect the value of contagion contact rates.



---

*Source: 5151674-2022-01-21.xml*