An efficient nonstandard computer method to solve a compartmental epidemiological model for COVID-19 with vaccination and population migration.

BACKGROUND AND
OBJECTIVE: In this manuscript, we consider a compartmental model to describe the dynamics of propagation of an infectious disease in a human population. The population considers the presence of susceptible, exposed, asymptomatic and symptomatic infected, quarantined, recovered and vaccinated individuals. In turn, the mathematical model considers various mechanisms of interaction between the sub-populations in addition to population migration.
METHODS: The steady-state solutions for the disease-free and endemic scenarios are calculated, and the local stability of the equilibium solutions is determined using linear analysis, Descartes' rule of signs and the Routh-Hurwitz criterion. We demonstrate rigorously the existence and uniqueness of non-negative solutions for the mathematical model, and we prove that the system has no periodic solutions using Dulac's criterion. To solve this system, a nonstandard finite-difference method is proposed.
RESULTS: As the main results, we show that the computer method presented in this work is uniquely solvable, and that it preserves the non-negativity of initial approximations. Moreover, the steady-state solutions of the continuous model are also constant solutions of the numerical scheme, and the stability properties of those solutions are likewise preserved in the discrete scenario. Furthermore, we establish the consistency of the scheme and, using a discrete form of Gronwall's inequality, we prove theoretically the stability and the convergence properties of the scheme. For convenience, a Matlab program of our method is provided in the appendix.
CONCLUSIONS: The computer method presented in this work is a nonstandard scheme with multiple dynamical and numerical properties. Most of those properties are thoroughly confirmed using computer simulations. Its easy implementation make this numerical approach a useful tool in the investigation on the propagation of infectious diseases. From the theoretical point of view, the present work is one of the few papers in which a nonstandard scheme is fully and rigorously analyzed not only for the dynamical properties, but also for consistently, stability and convergence.

Introduction

Epidemiology is considered a scientific discipline that studies the distribution, frequency, determinants, relationships, predictions, and control of factors related to health and disease in human populations [1]. Nowadays, epidemiology has a relevant place in many scientific areas, including the biomedical sciences, social sciences and even in the exact sciences [2]. In fact, it is worth pointing out that the study of diseases is an area which is as old as the birth of human writing. Indeed, the origins of the word “epidemiology” date back to ancient Greece, to some classical texts by Hippocrates of Kos, Aristotle and Galen [3]. Some of these scientists and philosophers were the first to use the terms “endemic” and “epidemic” in their works [4], though these concepts could have been used even before. However, epidemiology has witness a tremendous development since those times, being nowadays a useful discipline which encompasses various branches of human knowledge, even mathematical modeling and mathematical analysis [5]. These areas play an increasingly important role in the prediction and control of new pandemics like the coronavirus disease 2019 (SARS-CoV-2) or other diseases throughout human history [6].

It is important to recall that infectious diseases progress within populations due to both the behavior of the infectious agents and the population itself. Mathematical models which describe how an epidemic progresses are based on a set of assumptions and statistics that are used to establish suitable model parameters. In turn, these parameters completely determine the mechanics of propagation of the disease to a certain degree of reliability [7]. The mathematical models obtained in the way can be used then to predict which interventions to implement or avoid in order to control a disease, as well as patterns of growth and expansion that may result [8]. As expected, there is a vast amount of mathematical models which try to predict the evolution of a disease, and these models vary in complexity from simple deterministic models [9] to complex stochastic systems [10]. The former are usually based on differential or difference equations, while the latter employ usually stochastic equations. The approach chosen by epidemiologists depends on several variables including how much is known about the epidemiology of the disease, the purpose of the study, and the quantity and quality of data available [11].

Among the mathematical models used in mathematical epidemiology, compartment-based systems are a widely used technique for the quantitative and qualitative descriptions of the propagation of a disease [12]. This technique hinges mainly on describing the possible phases of interaction that a disease can have in a population [13]. It is worth pointing out that this type of models has been used to describe various diseases, and those systems are frequently based on the use of coupled ordinary differential equations [14]. Using this approach, several studies have been carried out to simulate the spreading of some diseases that have caused havoc in recent decades. For example, there are works which model and simulate the spread of Chikungunya disease [15], the control of measles in a human population [16], the epidemiology of diabetes mellitus with lifestyle and genetic factors [17], the epidemiology of sexually transmitted diseases [18], the modeling of tuberculosis disease in the Philippines [19], and the modeling of the coronavirus disease 2019 (COVID-19) pandemic [20], among other examples.

In the particular case of COVID-19, countries are currently working hard to fight this disease. To this day, this disease accounts for 5,732,354 deaths worldwide just 2 years after its first case [21]. Since then, many studies have been reported on the mathematical model of COVID-19, including some works using compartmental models to predict the effect of social distancing and vaccination as control measures [22], compartmental models for the COVID-19 pandemic with immunity loss [23], mathematical models for the calculation of COVID-19 lockdown efficiency [24] or the assessment of sensitivity and optimal economic evaluation with control intervention [25], a simple model without vaccination and migration [26], and even some compartmental models which employ various types of fractional-order operators in both space and time [27] among many examples available in the recent literature. In summary, various models have been proposed to describe the propagation of COVID-19 under various mathematical assumptions. It is worth mentioning that some of those works provide comparisons between various models and propose improvements in order to obtain more reliable paradigms. As an example, the authors of [28] carry out some detailed comparisons between various mathematical models and, after a careful analysis, they suggest that susceptible-exposed-infected-recovered-quarantined models are fundamental in order to capture the essential characteristics in the modeling of COVID-19.

The purpose of this work is to propose a general model that allows describing the spread of various diseases (including COVID-19) under general epidemiological assumptions. To that end, we will propose a compartmental system for an arbitrary human population. In particular, we will suppose that the population is separated into subpopulations of susceptible, exposed, symptomatic and asymptomatic infected, quarantined, recovered and vaccinated individuals. Various possible interactions between them will be taken into account, including the fact that recovered individuals may become susceptible. It is important to mention that the use of suitable model parameters will allow for the application of our mathematical model to particular diseases and epidemics. Our mathematical model will be based on the use of ordinary differential equations. We will determine the equilibrium points of this system along with their local and global stability properties, as well as the basic reproductive number. We will provide several simulations in this work, all of them obtained with a computer implementation of a nonstandard finite-difference method which is capable of preserving the most relevant analytical features of the solutions of the mathematical model. We must mention beforehand that the computational results will confirm the validity of our analytical properties. Finally, we will close this manuscript with a brief summary of the conclusions obtained in our study.

Methods

In this section, we deduce the mathematical model used to describe the propagation of a disease under suitable epidemiological assumptions. The epidemiological model will be analyzed to determine the equilibrium solutions and their stability properties. Among other analytical results presented in this section, we will derive the expression of the basic reproductive number using the next generation matrix approach [29].

To start with, we will consider a population of human individuals which are exposed to some contagious infection. Throughout, will represent the population size at the time , and we will suppose that the population is partitioned into the following seven compartments or subpopulations:

Susceptible individuals ().

Exposed individuals ().

Asymptomatic infected individuals ().

Symptomatic infected individuals ().

Quarantined individuals ().

Recovered/remove individuals ().

Vaccinated individuals ().

Obviously, the sizes of these subpopulations at time will be represented by , , , , , and , respectively. Under these assumptions, we have thatMoreover, to provide a more realistic epidemiological model, we consider in this work a constant migration into the population. More precisely, we will assume that a rate of people equal to , , and will migrate into the sub-populations of susceptible, exposed, asymptomatic and symptomatic.

Throughout this manuscript, all the parameters and variables in our mathematical model will take on non-negative real values. We will suppose that the population has natural birth and mortality rates which will be denoted by and , respectively. Susceptible individuals may become exposed if they have enough contact with exposed individuals at a rate of . On the other hand, susceptible individuals will be vaccinated a rate denoted by . Here, we will suppose that the vaccine is complete effective for all individuals, so it is appropriate to consider that vaccinated people will become susceptible at a rate equal to .

On the other hand, exposed individuals change compartment according to three possible options. The first one is to become quarantined, and we will assume that this will take place at a rate equal to . Alternatively, some exposed persons will become asymptomatic or symptomatic infected at rate equal to and , respectively. In turn, asymptomatic individuals may move to the recovered state at a rate given by . Individuals in the symptomatic compartment will become quarantined at a rate of . This may occur when the individuals present obvious symptoms of the disease. However, individuals can just move to the recovered state at a rate of or depending on whether then were quarantined or symptomatic. It is important to notice here that some quarantined and symptomatic individuals may die from the infectious disease, and we will employ and , respectively, to denote the rates at which these events occur. Finally, recovered individuals may become susceptible class with a rate equal to , under the assumption that the human body does not entirely create immunity to the disease. For convenience, Table 1 provides a summary of all the epidemiological parameters employed in this manuscript.

Table 1

Notations used in this work and their meaning.

Notations used in this manuscript and their meaning
Parameter	Description
Λ	Recruitment rate.
τ	Rate of transfer from vaccinated individuals to susceptible.
ω	Rate of transfer from susceptible individuals to vaccinated.
α	Contact rate between susceptible individuals and exposed individuals.
ζ	Rate of transfer of exposed individuals to quarantine.
ϵ	Rate of transfer of exposed individuals to symptomatic infected individuals.
δ	Rate of transfer of exposed individuals to asymptomatic infected individuals.
ι	Recovery rate of quarantine individuals.
υ	Mortality rate due to coronavirus in quarantine individuals.
κ	Rate of transfer of symptomatic infected individuals to quarantine.
ρ	Mortality rate due to coronavirus in symptomatic infected individuals.
θ	Recovery rate of transfer of symptomatic infected individuals.
σ	Rate of transfer of recovered individuals to susceptible.
μ	Natural mortality rate.
m_S	Rate of immigration of susceptible individuals.
m_E	Rate of immigration of exposed individuals.
m_I_A	Rate of immigration of asymptomatic infected individuals.
m_I_S	Rate of immigration of symptomatic infected individuals.

Figure 1 provides a flow chart which illustrates the epidemiological assumptions described above. Under these circumstances, the mathematical model describing the dynamics of propagation of the infectious disease is given by the following system of coupled nonlinear ordinary differential equations:The model will be complemented with initial conditions at the time . More precisely, we will assume that the initial compartment sizes will be provided by the non-negative numbers , , , , , and . Obviously, they will represent respectively the initial populations of susceptible, vaccinated, exposed, quarantined, asymptomatic infected, symptomatic infected and recovered.

Fig. 1

Flow chart describing graphically the dynamics of the compartmental epidemiological model proposed in this work.

It is important to notice that the mathematical model (2.2) has one disease-free equilibrium solution. To check this fact, let us assume a constant solution for the mathematical model in which . After some algebra, we readily check that the disease-free equilibrium is the point whose coordinates are given by

In order to calculate the basic reproductive number , we will employ the next generation matrix technique. Beforehand, recall that is the expected value of infection rate per time unit. Let us consider only those compartments of the mathematical model (2.2) which contribute to the dynamics of the infection, that is, let us consider the systemFollowing the approach in [29], we define the vectorsandTheir Jacobian matrices are, respectively,andA straightforward calculation shows thatwhere , and are real numbers, andAs a consequence, we obtain that the basic reproductive number is provided by the expression

Our next result summarizes the local stability analysis of the disease-free equilibrium.

The disease-free equilibrium of system(2.2)is locally asymptotically stable if.

Let represent the Jacobian matrix of associated to the system (2.2), and use to represent the matrix evaluated at the disease-free equilibrium solution. It is easy to check then that the Jacobian matrix is given bywhere Let be any complex number, and let be the identity matrix of order 7. If we let , then it is easy to check thatIn this expression for the matrix , we observe the following definitions for the components: Using properties of determinants, it is possible to check thatSetting the determinant equal to zero, solving for the unknown and rearranging terms algebraically, it follows that five of the eigenvalues of are The remaining eigenvalues satisfy the quadratic equationDescartes’ rule of signs imply that the number of negative roots for this polynomial is 0 or 2. However, the quadratic formula shows that the roots areIt follows that and are negative. Summarizing, notice that all the eigenvalues are negative if , in which case the disease-free equilibrium solution is locally asymptotically stable, as desired. □

Next, we proceed to calculate the endemic equilibrium solution. To that end, we assume a constant solution for the system (2.2), of the form , , , , , and , valid for all . Here, , , , , , and are non-negative constants. For the sake of convenience, we defineUnder these hypotheses, the mathematical model (2.2) reduces to the following system of algebraic equations:Proceeding algebraically, we may reach the identities: Moreover, after more tedious algebraic manipulations (or, equivalently, using symbolic software), it is possible to find out that We provide here the exact expressions for and only in view that they are relatively short. The expressions for the remaining coordinates of the endemic equilibrium point are actually too long to be written in this column. However, we must mention that all the coordinates are non-negative real numbers.

The endemic equilibrium point of system(2.2)is locally asymptotically stable if.

Notice that the Jacobian matrix evaluated at the endemic equilibrium point is given now bywhere However, the value of in the endemic case guarantees that . Using this fact and the properties of determinants, it is easy to check that the determinant of is given bywhere the expressions of the coefficients can be algebraically obtained, and are given in terms of the model parameters and the endemic equilibrium point. The exact expressions of these coefficients are long, and they were computed using symbolic algebra. We omit their expressions in view of the space available. The system is stable if the eigenvalues of the Jacobian matrix at the endemic point all have negative real parts. Using the Routh–Hurwitz criterion [30] and symbolic algebra, we obtain that the endemic equilibrium is stable if . □

In the next result, we will employ the gradient operator

The system(2.2)has no periodic solutions.

To establish this proposition, we will use the well known Dulac’s criterion. Let be the function defined component-wise for each by the expressionMoreover, letUsing differentiation, it is easy to check that We can check now thatThe conclusion follows now from Dulac’s criterion. □

Finally, we turn our attention to the problem on the existence and uniqueness of non-negative solutions of the mathematical model (2.2). To start with, it is obvious that solutions of (2.2) exist and are unique, for any set of initial conditions. This is a straightforward consequence of the fact that the model can be equivalently rewritten in the formwhere the function is given byfor each . Moreover, the function is given component-wise by , where each of the functions depends on , and is given by the right-hand side of the th differential equation in (2.2). The fact that is continuously differentiable assures the existence and uniqueness of continuous solutions for the mathematical model (2.2), for any set of initial conditions.

If the initial conditions,,,,,andare non-negative numbers, then the corresponding solution functions of the model(2.2)are likewise non-negative.

We proceed by contradiction. Suppose that some of the solutions take on negative values, and let be the greatest lower bound for which any of the solution functions is negative. Let represent the function for which this greatest lower bound occurs, and notice that . There are several cases whose proofs are entirely similar. We will only consider here the case in which . Observe then that and all the other functions at that time take on non-negative values. In particular, this implies that and . Using now that first equation of (2.2), it follows thatThus, there exists with the property that , for each . This contradicts the definition of , and we conclude that all the solution functions of the mathematical model (2.2) are non-negative for all times . □

Next, we would like to establish a bound for the growth of the population described by the epidemiological model (2.2). To that end, we add all the ordinary differential equations in (2.2) and simplify terms algebraically. It is easy to check thatAssuming that the initial population sizes of the compartments are non-negative, then the solution functions are likewise non-negative. As a consequence, the rate of change of increase of the total population is bounded from above by the non-negative constant . A straightforward integration yields then that, for each ,where . In view of this inequality, the following result is trivial.

Suppose that the initial conditions,,,,,andare non-negative numbers, and letbe a positive time period. Then the non-negative constantis a uniform bound for the solution functions of(2.2). □

As a consequence of the local stability properties, the boundedness of the solutions of the mathematical model (2.2) and the absence of periodic solutions, we conclude that the steady-state solutions are globally asymptotically stable.

Before closing this section, we provide a standard sensitivity analysis of the basic reproductive number with respect to the model parameters. To that end, for each parameter of the model, define the constantNotice then that andObserve that only , and are positive. We conclude that the basic reproductive number is sensitive only to the model parameters , and .

Results

In this section, we introduce a finite-difference scheme to approximate the solutions of (2.2). The methodology will be designed using the nonstandard approach popularized by Mickens in various of his seminal papers and monographs [31], [32], [33]. As the most important results, we will establish the main theoretical properties of our discretization, namely, the consistency, the stability and the convergence. Moreover, We will prove the capability of our scheme to preserve the positivity of the solutions, the constant solutions and their stability.

For the sake of convenience, agree that and , for each . We will approximate the solutions of our epidemiological model on a finite interval of time , where . Let , and fix a regular partition of the interval of the formfor each . For convenience, the associated partition norm will be represented by where, obviously, is a positive real number. We will use the lower-case symbols , , , , , , and to represent numerical approximations to the exact values of the functions , , , , , , and , respectively. Moreover, if is any of the lower-case symbols, then we will convey that , for each . Furthermore, we introduce the following linear discrete operator:It is well known that this operator provides a consistent approximation to the derivative if at the point , with consistency order equal to one in time. Alternatively, it also yields a first-order consistent approximation of the derivative of with respect to at the time .

Using this nomenclature, the finite-difference scheme employed to approximate the solutions of the system (2.2) at time is given by the algebraic nonlinear system of equationsObviously, this is a nonstandard discretization in the sense that the approximations provided for some terms in the scheme are provided in a non-local manner. The numerical model is a two-step system which will be theoretically analyzed in this section. To that end, it is important to notice that the discrete model (3.3) can be alternatively expressed in explicit form. After some algebraic manipulations, the finite-difference scheme can be equivalently rewritten as

From this discussion, it is obvious that the discrete model (3.3) is a semi-explicit algebraic system. We just need to point out that is given in terms of as the first equation of (3.4) shows. However, this shortcoming can be saved calculating firstly from the third equation, and the obtaining from the first identity of (3.4). Moreover, the following theoretical result is also straightforward. The reader will notice that this is the discrete version of Theorem 4.

If the initial conditions,,,,,andare non-negative, then the discrete system(3.3)has a unique solution, and all the solution functions are non-negative.

If , then , , , , , and are non-negative numbers by hypothesis. Now, suppose that the conclusion of this result is true for some . Under these circumstances, the right-hand sides of the identities in (3.4) are non-negative. As a consequence, the approximations at time are also non-negative, and the conclusion of this theorem follows by induction. □

Our next step is to obtain a discrete form of the inequality (2.72). To that end, let us suppose that the initial conditions are all non-negative. As a consequence of the previous theorem, the numerical solutions are likewise non-negative. Add together the equations in the discrete system (3.3) and simplify terms. It is easy check that if , thenLet , and take the sum on both sides of this equation for between 0 and . Using the formula for telescoping sums, simplifying algebraically, rearranging terms, recalling that the solutions of the discrete model (3.3) are non-negative and using the fact that , we obtain the following upper bound for the total population at the time :Observe that the continuous estimate (2.72) is recovered from this last inequality when we let . Moreover, we have the following discrete version of Theorem 5.

Suppose that the initial conditions,,,,,andare non-negative. Then the non-negative numberis a uniform bound for the solutions of model(3.3). □

The following theorem establishes that the disease-free and the endemic equilibrium solutions are also constant solutions of the numerical model (3.3). Moreover, their stability properties are also preserved in the discrete scenario.

The pointsandare constant solutions of the numerical model(3.3). Moreover, the following hold:

The pointis locally asymptotically stable if, and unstable if.

The pointis locally asymptotically stable if, and unstable if.

The points are constant solutions of (3.3) follows after a simple substitution in that system. On the other hand, the local stability properties of the numerical model are satisfied in view that the Jacobian matrix of the discrete system (3.3) is the same as that of the continuous model (2.2). □

Next we establish the numerical properties of our finite-difference scheme. More precisely, we will prove that the numerical scheme (3.3) is consistent, stable and convergent. To that end, let us define the differential operatorsObviously, these continuous operators are defined for each . For each and being any of the solution functions of (2.2), agree that . On the other hand, define also the difference operators

Let us define now for each . In the following, we will use and to denote, respectively, the -norm and the infinity norm in . More precisely, if , then Moreover, we introduce the norm

Using this nomenclature, we will prove firstly the consistency of the finite-difference method (3.3).

Ifare of class, then there exists a constantwhich is independent of, such that.

Let . Using the regularity of the function , it follows that this function is bounded on by some non-negative constant . Moreover, since is of class , Taylor’s theorem readily guarantees that there exist constants which are independent of , with the property that for each . As a consequence of these inequalities, algebraic simplifications and the triangle inequality, it follows thatfor each . Here,which is a non-negative constant that is independent of . It follows then that . In similar fashion, it is possible to show that there exist constants which are independent of for , such that the inequality is satisfied. Now, if we define the non-negative constantthen the constant is independent of and it satisfies the conclusion of this theorem, as desired. □

Next, we turn our attention to the stability and convergence properties of the finite-difference scheme (3.3). The following discrete form of Gronwall’s inequality will be needed.

Pen-Yu [34]

Letandbe finite sequences of nonnegative mesh functions, and suppose that there existssuch thatThenfor each. □

To establish the stability property, we will consider two sets of non-negative initial conditions for the finite-difference scheme, which we will denote respectively by According to Theorem 6, the discrete model (3.3) yields non-negative solutions for each of these solutions. These solutions will be denoted respectively by and , where Moreover, we will agree that , for each and . This nomenclature will be used in the following theorem.

Letandbe non-negative initial conditions for the model(3.3), and suppose thatandare the respective solutions. Ifis sufficiently small, then there is a constantsuch that, for each. Here,

For the sake of convenience, we will let , for each . It is obvious thatis satisfied for each . Moreover, after simplification and some additional algebraic steps, it is possible to show the sequence satisfies the system of algebraic equationsLet , and assume that is the uniform bound for the solution functions of (3.3). Take absolute value on both sides of the first equation of the discrete system (3.27), groups similar terms algebraically and use the triangle inequality to obtainfor each . Here, . Multiplying both sides of this inequality by , summing on both sides for all from 0 to , using the formula for telescoping sums and rearranging terms, we obtain thatIn turn, if we letthenfor each . In similar fashion, we may use the remaining equations in (3.27) to show that there exist non-negative constants and for each , such that and

Let be sufficiently small so that , for each , and let satisfy , for each . Adding the inequalities (3.31)–(3.37) and lettingwe obtain thatFinally, we use Lemma 1 to establish that , where , for each . The conclusion of this theorem readily follows from this fact. □

In terms of the nomenclature employed in the proof of the previous theorem, observe that the conclusion can be rewritten as , for each .

Our final theoretical result summarizes the convergence property of the finite-difference method (3.3). We omit the proof in view that it is similar to that of Theorem 10. We just need to point out that is the difference between the exact solution and the numerical approximation , for each and . The consistency property of the computer method is also required to bound the local truncation error, along with the discrete form of Gronwall’s inequality.

Suppose that the solutions of problem(2.2)are of class. For sufficiently small values of, the solutions of the discrete model(2.2)converge in the-norm to the exact solution with order of convergence equal to. □

Before closing this section, we present computer simulations which confirm the validity of some of the analytical results derived in this work. Our simulations have been carried out with the Matlab code provided in Appendix A. It is worth pointing out that the computational implementation is relatively simple, which is yet another important advantage of the approach introduced in the present manuscript. It is worth pointing out that the parameter values will be those in Table 2, and that some of those values wer taken from [35].

Table 2

Values of the parameters used in the various computational experiments presented in this manuscript.

Parameter	Value
α	0.01
δ	1.6728×10^−5
ϵ	0.0101
ζ	0.02798
η	0.04478
θ	0.0101
ι	0.0045
κ	0.0368
Λ	0.06
μ	0.0106
ρ	0.004
σ	0.0668
τ	0.0002
υ	3.2084×10^−4
ω	0.0032
m_S	0
m_E	0
m_I_A	0
m_I_S	0

In our first example, we will confirm the local stability properties of the disease-free equilibrium solution. To that end, we will employ the parameter values in Table 2 along with the data set 1 from Table 3. Under these circumstances, Figure 2 shows the dynamics of the solution for (a) the susceptible population and (b) the vaccinated population with respect to time, over the time period . The results confirm the stability of the disease-free equilibrium solution. It is worth pointing out that the value of the basic reproductive number is equal to 171.12. Moreover, the dashed line represents the value of the equilibrium point. Obviously, the solutions tend to reach those values as tends to infinity. In turn, (c) and (d) show, respectively, the susceptible population and the vaccinated individuals as functions of time, for . We employed data set 2 from Table 3 in this case. Again, these sub-populations converge asymptotically to the steady-state solutions. So, whether the initial conditions are close to or far from the equilibrium solutions of the system (2.2), the numerical solutions converge to these values as time increases. This is in agreement with the analytical results. Moreover, the simulations show that the numerical method preserves the equilibria and their stability, as expected from the theoretical analysis. □

Table 3

Initial conditions used in the numerical experiments of this manuscript.

	Data
Parameter	Set 1	Set 2	Set 3
S^0	S_0−0.3	S_0+30	20
V^0	V_0+0.3	0	0
E^0	0	0	1
IA0	0	0	0
Q^0	0	0	0
IS0	0	0	0
R^0	0	0	0

Fig. 2

Graphs of the temporal behavior of the sub-population of susceptible (first column) and the sub-population of vaccinated (second column) in the mathematical model (2.2). We employed the parameter values in Table 2, along with the initial conditions given by data set 1 (first row) and data set 2 (second row) in Table 3. The dashed lines represent the theoretical steady-state solutions for the disease-free scenario.

In this example, we consider the endemic case and show once more that the finite-difference scheme is capable of preserving the steady-state solutions and their stabilities. Moreover, we provide computational proof that the endemic equilibrium is globally asymptotically stable as proved in the previous section. To that end, consider the parameter values given in Table 2, along with the initial conditions under data set 3 of Table 3. The results are provided in Figure 3 as time-dependent graphs of (a) susceptible, (b) vaccinated, (c) exposed and (d) asymptomatic infected, and in Figure 4 as graphs of (a) quarantined, (b) symptomatic infected, (c) recovered and (d) total population, for . For convenience, the theoretical endemic equilibrium values are plotted as dashed lines. The results show that the solutions tend to their equilibrium values as time increases. This is in obvious agreement with the theoretical results derived in this work. □

Fig. 3

Graphs of the temporal behavior of the sub-populations of (a) susceptible, (b) vaccinated, (c) exposed and (d) asymptomatic infected individuals in a population modeled by (2.2). We employed the parameter values in Table 2, along with the initial conditions given by data set 3 3. The dashed lines represent the theoretical steady-state solutions for the endemic scenario.

Fig. 4

Graphs of the temporal behavior of the sub-populations of (a) quarantined, (b) symptomatic infected, (c) recovered and (d) total population of individuals in a system modeled by (2.2). We employed the parameter values in Table 2, along with the initial conditions given by data set 3 3. The dashed lines represent the theoretical steady-state solutions for the endemic scenario.

Conclusions

In this work, we investigated both analytically and numerically a compartmental epidemiological model which describes the propagation of a disease among a human population. The model is intended to describe the propagation of COVID-19, but it can be used to any other disease which satisfies the epidemiological hypotheses used in this work. Among the distinctive features of the model, we considered various compartments: susceptible, exposed, asymptomatic and symptomatic infected, quarantined, recovered and vaccinated individuals. We supposed also that population migration is possible in the mathematical model. Analytically, we obtained the steady-state solutions of the model, and determined conditions for their local stability. The basic reproductive number was determined using the next generation matrix, and we established the existence and uniqueness of non-negative solutions. Also, we provided an upper bound for the solutions functions, and the analysis of parametric sensitivity was theoretically carried out.

As one of the most important results of our study, we proposed a finite-difference method to approximate the solutions of our mathematical model. This computational technique was designed using the nonstandard approach proposed by R. E. Mickens. An explicit form of the scheme was provide, and we established the existence and uniqueness of non-negative solutions for this mathematical model. We showed that the computational scheme has the same steady-state solutions as the continuous model and, moreover, the stability properties are also preserved by our discretization. We proved that the discrete model is a consistent discretization of the epidemiological model, and the conditional stability and convergence properties were derived using a discrete form of Gronwall’s inequality. Here, it is worth pointing out that many nonstandard techniques are usually presented in the literature without providing these numerical properties. However, we established them mathematically in the present manuscript. Computationally, we obtained various simulations to illustrate the performance of our scheme. The results showed that the method identifies correctly the steady-state solutions and, moreover, it is also able to reproduce the stability properties of the continuous model. Our simulations show additionally that the scheme is capable or preserving the non-negativity and the boundedness of the solutions, in agreement with out theoretical results.

It is important to point out that the discretization proposed in this work is first order accurate in the temporal variable. As one of the reviewers pointed out, this numerical accuracy may not be satisfactory in the practice, in particular when dynamical simulations are performed. Other approaches may have the advantage of providing an accuracy of higher order, like the family of Runge–Kutta methods for systems of ordinary differential equations, which is a family of stable and convergent techniques. However, those approaches may not be able to preserve the positivity and boundedness of the solutions, or may not be able to preserve the equilibria and their stability. Nevertheless, in the case that they can preserve those features, the present methodology has its simplicity as one of the advantages. As the appendix shows, the present methodology is relatively easy to implement even for a scientist with little knowledge in computer programming.

On the other hand, as one of the anonymous reviewers of this manuscript pointed out, it is important to mention that there exist various reports available in the literature of third-order methods for time-dependent nonlinear partial differential equations in which the convergence and the stability have been analyzed. For example, there are reports on fully discrete Fourier collocation spectral methods for the 3-D viscous Burgers equation [36], high-order multi-step numerical schemes for two-dimensional incompressible Navier–Stokes equations [37], high-order exponential time-differencing numerical schemes for no-slope-selection epitaxial thin-film models with energy stability [38], third-order BDF energy-stable linear schemes for the no-slope-selection thin film model [39] and BDF-type energy-stable schemes for the Cahn–Hilliard equation [40].

Finally, as the same reviewer pointed out, there are various reports in which some logarithmic energy potential has been introduced for reaction-diffusion equations or other related gradient flows. In such way, the preservation of the positivity of the solutions has been ensured. As examples, we can mention some numerical works for the Poisson–Nernst–Planck system [41], a ternary Cahn–Hilliard system with the singular interfacial parameters [42], the three-component Cahn–Hilliard-type model for macromolecular microsphere composite hydrogels [43], the binary fluid-surfactant system [44], a liquid thin-film coarsening model [45], the Poisson–Nernst–Planck–Cahn–Hilliard equations with steric interactions [46], the Cahn–Hilliard equation with variable interfacial parameters [47], the Cahn–Hilliard equation with a Flory–Huggins–Degennes energy [48], he Cahn-Hilliard equation with logarithmic potential [49] and a reaction-diffusion system with detailed balance [50], [51]. The authors of the present manuscript are not aware whether an entropy could be introduced in the epidemiological mathematical model (2.2), in order to guarantee the preservation of the positivity of the numerical solutions.

Funding

One of the authors (J.E.M.-D.) was funded by the National Council of Science and Technology of Mexico through grant A1-S-45928.

Data statement

The data that support the findings of this study are available from the corresponding author, J.E.M.-D., upon reasonable request.

Conflict of interest statement

The authors declare no potential conflict of interest.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.