A Discrete-Time Compartmental Epidemiological Model for COVID-19 with a Case Study for Portugal

In [Ecological Complexity 44 (2020) Art. 100885, DOI: 10.1016/j.ecocom.2020.100885] a continuous-time compartmental mathematical model for the spread of the Coronavirus disease 2019 (COVID-19) is presented with Portugal as case study, from 2 March to 4 May 2020, and the local stability of the Disease Free Equilibrium (DFE) is analysed. Here, we propose an analogous discrete-time model and, using a suitable Lyapunov function, we prove the global stability of the DFE point. Using COVID-19 real data, we show, through numerical simulations, the consistence of the obtained theoretical results.


Introduction
The coronavirus belong to the family of Coronaviridae. It is a virus that causes infection to humans, other mammals and birds. The infection usually affects the respiratory system, and the symptoms may vary from simple colds to pneumonia. To date, eight coronaviruses are known. The new coronavirus SARS-CoV-2 originates COVID-19, and was identified for the first time in December 2019 in Wuhan, China. Two other coronaviruses have caused outbreaks: SARS-CoV, in 2002-2003, and MERS-CoV in 2012 [1]. The COVID-19 pandemic in an ongoing pandemic, probably the most significant one in human history. The effects are so severe that it has been necessary to use quarantining to control the spread. It has forced humans to adapt to this new situation [2].
The use of mathematical compartmental models to study the spread and consequences of infectious diseases have been used successfully for a long time. The techniques to study compartmental models are vast: continuous models, fractional models, and discrete models, among others [2]. The complexity of the models is determined by the effects and consequences of the disease. Several models have been proposed for COVID-19 spread; see [2][3][4][5], among others. In this work, we are interested in discretizing the model presented in [3], using the Mickens method [6,7], which was applied successfully in several different contexts [8].
In Section 2, we recall the SAIQH (susceptible-asymptomatic-infectious-quarantinedhospitalized) continuous-time mathematical model, the equilibrium points and the stability results of [3]. Section 3 is dedicated to the original results of our work: the new discrete-time model, the proof of its well-posedness, the computation of the equilibrium points, and the proof of the stability of the disease-free equilibrium point. We end Section 3 with numerical simulations showing the consistency of our results. Finally, Section 4 is dedicated to discussion of the obtained results and some possible future work directions. arXiv:2111.11860v1 [math.NA] 23 Nov 2021

Preliminaries
This section is dedicated to the presentation of the continuous-time model [3] and its results, which establish the stability of the equilibrium points. For more information, see [3].
The total living population under study at time t ≥ 0 is denoted by N(t) as follows: where S(t) represents the susceptible individuals; A(t) the infected individuals without symptoms or with mild ones; I(t) the infected individuals; Q(t) the individuals in quarantine, that is, in isolation at home; H(t) the hospitalized individuals; and, finally, H(t) the hospitalized individuals in intensive care units. There is another class that is also considered, D(t), that gives the cumulative number of deaths due to COVID-19 for all t ≥ 0. Regarding the parameters, all of them are positive. The recruitment rate into the susceptible class is Λ; all individuals of all classes are subject to the natural death rate µ along all time t ≥ 0 under study. Susceptible individuals may become infected with COVID-19 at the following rate: where β is the human-to-human transmission rate per unit of time (day) and l A and l H quantify the relative transmissibility of asymptomatic individuals and hospitalized individuals, respectively. The class H does not enter in λ because the number of health care workers that become infected by SARS-CoV-2 in intensive care units is very low and can be neglected. A fraction p ∈ [0, 1] of the susceptible population is in quarantine at home, at rate φ. Consequently, only a fraction 1 − p ∈ [0, 1] of susceptible individuals are assumed to be able to become infected.
Since there is uncertainty about long-immunity after recovery, it is assumed that individuals of class Q will become susceptible again at a rate ω. It is also considered that only a fraction m ∈ [0, 1] of the quarantined individuals move from class Q to S. It means that (m × 100)% of the quarantined individuals return to class S at the end of 1 ω days. These assumptions are justified by the state of calamity that was immediately decreed by the government of Portugal to address the state of emergency, which was fully respected by the Portuguese population. After 1 ν days of infection, only a fraction q ∈ [0, 1] of infected individuals without (or with mild) symptoms have severe symptoms. Thus, (q × 100)% of individuals of compartment A move to I at rate ν. A fraction f 1 ∈ [0, 1] of infected individuals with severe symptoms are treated at home and the other fraction (1 − f 1 ) ∈ [0, 1] are hospitalized, both at rate δ 1 . The model considers the three following scenarios for hospitalized individuals:

1.
A fraction f 2 ∈ [0, 1] of individuals in class H can evolve to a state of severe health status, needing an invasive intervention, such as artificial respiration, so they need to move to intensive care, at rate δ 2 ; 2.
A fraction f 3 ∈ [0, 1] of individuals in class H die due to COVID-19, the disease related death rate associated with hospitalized individuals being α 1 ; 3.
Regarding the hospitalized individuals in intensive care units, the model considers two possibilities as follows:
A fraction κ ∈ [0, 1] of individuals in class H die due to COVID-19, the disease-related death rate associated with hospitalized individuals in intensive care units being α 2 .
Compiling all the previous assumptions, one has the following mathematical model: It can be seen in [3] that (3) is equivalent to the following: Table 1 presents all parameters and initial conditions of mathematical model (4).
In [3], it is shown that the biologically feasible region is given by the following: which is positively invariant for (4) for all non-negative initial conditions. To simplify the expressions in the computations, the following notation is used: Model (4) has the following reproduction number: and two equilibrium points, the disease free equilibrium (DFE) point and the endemic equilibrium (EE) point that is given by the following: Assuming that the transmission rate is strictly positive, we have the following: Using (8) in (9), it can be seen that the endemic equilibrium satisfies the following: Regarding the stability of the equilibrium points, the following results hold.
). The disease free equilibrium point Σ 0 (7) is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1, where R 0 is the basic reproduction number (6).
The model (4) has a unique endemic equilibrium point whenever R 0 > 1.

Results
In this section, we begin by presenting our proposal for a COVID-19 discrete-time model. After that, we show the well-posedness of the model, that is, we prove that the solutions of the model are positive and bounded and that the equilibrium points coincide with those of the continuous-time model presented in Section 2. We finalize this section by proving the global stability of the disease-free equilibrium point, using a suitable Lyapunov function and with the presentation of some numerical simulations, using real data, showing the consistency of our theoretical results.

The Discrete-Time Model
One of the important features of the discrete-time epidemic models obtained by the Mickens method is that they present the same features of the corresponding original continuoustime models. Our nonstandard finite discrete difference (NSFD) scheme for solving (4) is a numerical dynamically consistent method based on [9]. Let us define the time instants t n = nh with n integer, the step size as h = t n+1 − t n , and (S n , A n , I n , Q n , H n , H n ) as the approximated values of the following: (S(nh), A(nh), I(nh), Q(nh), H(nh), H(nh)).
Discretizing system (4) using the NSFD scheme, we obtain the following: 6 of 15 where the denominator function is ψ(h) = exp(µh)−1 µ [10]. Throughout this work, for brevity, we write ψ(h) = ψ. Remark 1. In the continuous model, the reproduction number (6) can be rewritten as the following: Let us consider the region Ω = (S, A, I, Q, H, H) ∈ (R + 0 ) 6 : 0 < N ≤ Λ µ . Our next lemma shows that the feasible region is equal to the continuous case.  (11) is linear in (S n , A n , I n , Q n , H n , H n ), we can rewrite it as the following:

Proof. Since model
Since all parameters of model (13) are positive and the initial conditions are also positive, then, by induction, S n ≥ 0, A n ≥ 0, I n ≥ 0, Q n ≥ 0, H n ≥ 0 and H n ≥ 0 for all n ∈ N. Regarding the boundedness of the solutions, the total population N n = S n + A n + I n + Q n + H n + H n . Adding all equations of (13), we obtain the following: By Lemma 2.2 of Shi and Dong [11], then N n ≤ Λ µ for all n ∈ N. Thus, Ω is the biologically feasible region.

Global Stability
We now prove the global stability of (14).
Theorem 1. For the discretized system, the DFE point E 0 is globally stable if R 0 < 1. If R 0 > 1, then the DFE is unstable.
Proof. Let us define the discrete Lyapunov function L n as the following: L n (S n , A n , I n , Q n , H n , H n ) = 1 Hence, L n (S n , A n , I n , Q n , H n , H n ) ≥ 0 for all ξ n ≥ 0, ξ ∈ S, A, I, Q, H, H . Moreover, L n (S n , A n , I n , Q n , H n , H n ) = 0 if and only if (S n , A n , I n , Q n , H n , H n ) = E 0 . Computing ∆L n = L n+1 − L n , we have the following: From equations of system (11), we obtain the following: and, at the disease-free equilibrium point E 0 , we have the following relations: Substituting in ∆L n , and simplifying, we obtain the following: Then, From the equations of (13), it can be seen that and Since is equal to the following: Gathering the information from (21) and (23), (20) becomes the following: Once more, from system (13) it can seen that −µ(A n+1 + I n+1 ) < −a 0 A n , so Hence, if R 0 < 1, then ∆L n ≤ 0 for all n ≥ 0, that is, L n is a monotone decreasing sequence. We have L n ≥ 0. Then there is a limit for lim n→∞ L n ≥ 0. Therefore, lim n→∞ ∆L n = 0 implies that lim n→∞ S n = S 0 , lim n→∞ Q n = Q 0 , lim n→∞ A n = lim n→∞ I n = lim n→∞ H n = lim n→∞ H n = 0. So, if R 0 < 1, then E 0 is globally asymptotically stable.

Numerical Simulations
In this section, we show, numerically, that our discretized model describes well the transmission dynamics of COVID-19 in Portugal, from 2 March to 4 May 2020. Our data were obtained from the daily reports from DGS, available in [12], and since we consider the spread from 2 March, the initial time t = 0 corresponds to 2 March 2020.
The initial values are the same as those used in [3]. Regarding the parameters, using [13], we can obtain the values from 2019, so the parameters were updated with respect to [3], but the difference is very small. All the values used in our simulations are presented in Table 2. Table 2. Parameter values and initial conditions of (4).  [12] Using the values of Table 2, in Figure 1, we compare the number of infected individuals predicted by our discrete-time model, the ones predicted by the continuous model of [3], and the real data.   Table 2, R 0 ≈ 0.95, the number of infected and asymptomatic individuals, as well as the ones in the hospital and intensive care units, tend to disappear. The susceptible individuals and those in quarantine tend to the equilibrium point E 0 as time increases.

Real Data Inf-Model
Comparing the graphics from Figures 2-7, we can see that the asymptotic behavior is similar. However, the discrete-time model predicts, in some instances, smaller values. All the computations were done using Mathematica, version 12.1.

Discussion
The COVID-19 pandemic has been a shocking experience to every human being, and no one expected the consequences felt worldwide. From this, we can also state that the number of variables and parameters to be taken into account in a COVID-19 model are not small. Indeed, if there is such thing as a good model that explains COVID-19, surely it is not a simple model. Moreover, most of the models that explain epidemiological, ecological, and economical phenomena are not linear. Therefore, given the tools available, it is impossible to present an exact solution to the problem. That is one of the reasons for the development of nonstandard finite difference methods. Indeed, such methods enable us to construct discrete models from continuous ones, allowing one to present numerical solutions.
Nowadays, there are several nonstandard finite schemes. We have used here the one developed by Mickens because it has shown to be dynamically consistent. Precisely, we discretized the continuous-time model presented and analyzed in [3]. With our discrete-time model, we analyzed the evolution of the COVID-19 pandemic in Portugal, from its beginning up to 4 May 2020. This choice allowed us to compare our results with those of [3]. We obtained the equilibrium points of the proposed discrete-time model, showing that they coincide with the ones from the continuous model. From the qualitative analysis point of view, we believe that it is important to prove the stability of the equilibrium points if the model describes real data, which is our case here. For this reason, Figures 2-7 show the simulation results made over 2000 days: our intention is to show the stability of the model. Regarding the stability, it should be noted that, in contrast with [3], which only proves local stability, here, we have proved the global stability of the disease-free equilibrium (DFE) and attempted to do the same with the endemic equilibrium (EE) point. However, regarding the EE, we were not able to prove the global stability analytically, and the question remains open. Figure 1 compares the predictions of the continuous model with those of the discrete-time model. One concludes that the discrete-time model fits a little bit better to the real data. From Figures 2-7, we present the convergence to the disease-free equilibrium point, using both continuous and discrete models. The asymptotic behavior is similar; the only difference worth mentioning is that the discrete-time model, in some small interval of time, predicts smaller values than the continuous one.
In this work, we restricted ourselves to the development of the pandemic in Portugal. We are aware that new models are emerging and, in a future work, we intend to analyze a model that fits the data of more than one country and region, possibly globally. This is under investigation and will be addressed elsewhere. Another line of research concerns models that take vaccination into account [18]. Here, vaccination is not considered because our model was created to explain the development of COVID-19 at the beginning of the pandemic. However, regarding vaccination, individuals do not become immune, so we think that we need to wait a little longer to see the development. It remains open the question of how to prove global stability for the endemic equilibrium.