key: cord-0853377-upf3vqe8
authors: Rafiq, Muhammad; Macías-Díaz, J.E.; Raza, Ali; Ahmed, Nauman
title: Design of a nonlinear model for the propagation of COVID-19 and its efficient nonstandard computational implementation
date: 2020-09-22
journal: Appl Math Model
DOI: 10.1016/j.apm.2020.08.082
sha: 3d05308548716418e26d09a1ef1dba52f8c26600
doc_id: 853377
cord_uid: upf3vqe8

In this manuscript, we develop a mathematical model to describe the spreading of an epidemic disease in a human population. The emphasis in this work will be on the study of the propagation of the coronavirus disease (COVID-19). Various epidemiologically relevant assumptions will be imposed upon the problem, and a coupled system of first-order ordinary differential equations will be obtained. The model adopts the form of a nonlinear susceptible-exposed-infected-quarantined-recovered system, and we investigate it both analytically and numerically. Analytically, we obtain the equilibrium points in the presence and absence of the coronavirus. We also calculate the reproduction number and provide conditions that guarantee the local and global asymptotic stability of the equilibria. To that end, various tools from analysis will be employed, including Volterra-type Lyapunov functions, LaSalle’s invariance principle and the Routh–Hurwitz criterion. To simulate computationally the dynamics of propagation of the disease, we propose a nonstandard finite-difference scheme to approximate the solutions of the mathematical model. A thorough analysis of the discrete model is provided in this work, including the consistency and the stability analyses, along with the capability of the discrete model to preserve the equilibria of the continuous system. Among other interesting results, our numerical simulations confirm the stability properties of the equilibrium points.

Coronavirus disease 2019 (COVID-19) is a viral disease that was identified toward the end of the year 2019 in China, and which became a pandemic in the first quarter of the year 2020 [1] . After its identification in 2019, COVID-19 has been a source of active research from various scientific points of view, mainly due to the mortality rate derived from this disease and the inherent health complications. Indeed, to this day, the COVID-19 dashboard by the Center for Systems Science and Engineering at Johns Hopkins University reports a total of 5,439,559 cases in 188 countries around the World, along with 2,185,696 recovered individuals and 345,589 deaths worldwide [2] . On the other hand, the usual symptoms of this disease include fever, cough, fatigue and shortness of breath [3] , but it may develop to acute respiratory distress syndrome, multi-organ failure, septic shock and blood clots, among other conditions [4, 5] . In some cases, these symptoms evolve rapidly resulting ultimately in a painful death. Needless to mention that there are other severe consequences on the propagation of this disease worldwide, including economic and financial outcomes which are difficult to assess at this stage of the pandemic. Some reports already have started to investigate these potential consequences at various geographic levels [6, 7] , though the full extent of the implications of the global spread of COVID-19 are still under scientific scrutiny.

As many other problems associated to the study of COVID-19, the reliable prediction of its dynamics is an important topic of investigation from the epidemiological point of view. In general, the dynamics of epidemics is a highly transited topic of research in mathematics and simulation. There are various families of models depending on the assumptions on the mechanisms of propagation. For example, there are susceptible-infected-susceptible (SIS) systems which have been used to model nonlinear incidence rates and double epidemic hypotheses [8] , susceptible-infectedrecovered (SIR) models which are employed to consider the awareness of the presence of a disease [9] , susceptibleexposed-infected-recovered (SEIR) models that predict the propagation of epidemics with non-local reaction functions [10] , susceptible-exposed-infected-quarantined-recovered (SEIQR) systems used to account for adjusted incidence and imperfect vaccinations [11] , among other examples. It is worth pointing out that most of the deterministic models available in the literature are based on the use of ordinary differential equations, which means that a constant diffusion among the population is assumed. However, there are approaches which hinge on the use of systems of partial differential equations in order to account for imperfect mixing modelled by diffusion [12, 13, 14, 15] . In the present work, we focus our attention to the spreading of COVID-19 in populations where the diffusion is constant.

It is important to realize that the most recent progresses in mathematical epidemiology give account of a vast variety of realistic mathematical models to describe the propagation of diseases [16] . Indeed, in addition to the usual models based on the use of ordinary or partial differential equations, there are reports on individual-based systems to describe the spreading of infectious diseases in clustered complex networks [17] , agent-based systems to investigate the spatio-temporal dynamics and synchrony of influenza epidemics [18] , network-based models to provide control strategies in SIS systems [19] , stochastic models that consider ratio-dependent incidence rates [20] , statistical models for the problem of estimation in some epidemic scenarios [21] , among various other approaches. The analysis of those models has been a fruitful topic of research in the literature. Moreover, the reliable computational simulation of those systems has also become a cornerstone in order to propose pertinent applications [22] . In the case of models based on differential equations, consistency, stability and convergence of the numerical methodologies have been criteria for the reliability of the discretizations. However, in light that the variables under study are populations, emphasis has been done on the development of finite-difference schemes that preserve the positivity and/or boundedness of the solutions of epidemic problems [23, 24, 25] , as well as other approaches for more complicated scenarios [26] .

In this manuscript, we will consider a SEIQR model for the propagation of COVID-19. The system will consist of a nonlinear set of coupled ordinary differential equations in which various real parameters are included. The use of such model is based on the current understanding of the mechanisms of propagation of COVID-19. In particular, the consideration of the subpopulations of susceptible, exposed, infected, quarantined and recover individuals is strongly supported by both the effects of COVID-19 in human health and the measures adopted by the governments of all countries affected [27] . Various assumptions will be imposed on the problem considered in this work, and we will obtain a system of ordinary differential equations from those hypotheses. We will provide a stability analysis of the equilibrium points of the system in the presence and absence of COVID-19. The nature of the local stability of those equilibria will be thoroughly investigated with respect to the basic reproduction number. As a means to verify our theoretical results, we will propose a two-level finite difference discretization of the SEIQR model. Our approach will hinge on a non-local discretization of the continuous problem, and we will show that the discrete model has various numerical and structural properties which make it a desirable tool in the simulation of the SEIQR system. For instance, we will prove the linear stability of the scheme, along with its capability to preserve the positivity of solutions. This feature is of the utmost importance in view that the variables under investigation are population sizes [28, 29] . Some simulations will be provided to assess the validity of the theoretical results derived in this work.

The present work is organized as follows. In Section 3, we will introduce the mathematical model under investigation. As we pointed out before, the system will be a general nonlinear SEIQR model with constant spatial diffusion. The assumptions on the derivation of the model will be provided therein, and the model will include various physically relevant constants to make it as realistic as possible. In that section, we will obtain the equilibrium points both in the presence and the absence of COVID-19. We will obtain the reproduction number in that stage, and various results on the local stability of the equilibria will be rigorously established. The numerical model employed to solve the SEIQR model is presented in Section 3. The scheme is a nonstandard model in the sense of R. E. Mickens [30] . We will show therein that the scheme is a stable technique which is able to preserve the positivity of the solutions. In other words, our computational methodology is a structure-preserving technique or, equivalently, a dynamically consistent technique. As the end of that section. we provide some simulations that illustrate the main analytical features of the numerical methodology, and we close this manuscript with a section of concluding remarks. 

For the remainder of this work, we will investigate the dynamics of a human population which is divided into various non-overlapping subpopulations. The population will be represented by P : [0, ∞) → R, and it will be a function of the time t ≥ 0. In turn, each of the five subpopulations will be denoted by the nonnegative differentiable functions S , E, I, Q, R : [0, ∞) → R, which will represent, respectively, the total number of individuals who are susceptible to be infected by COVID-19, the number of individuals who have been exposed to this disease, the number of individuals who are infected, the number of persons who are in quarantine and the number of recovered individuals.

In this work, we will assume that the population of susceptible individuals is increased constantly as a recruitment rate λ. At the same time, we will consider a natural death rate µ of susceptible individuals. Moreover, the population of susceptible people will decrease at constant rates when they interact with exposed and with infected individuals. Those constant rates are β 1 and β 2 , which represent the contact rate of susceptible individuals with infected and exposed individuals, respectively. Meanwhile, the population of exposed individuals increases at a rate which is proportional to the interactions between the susceptible population and both the infected and exposed populations. The respective constants of proportionality are equal to β 1 and β 2 . On the other hand, the rate of change of exposed also decreases as a factor of the quarantined rate of exposed individuals q 1 , the recovery rate of exposed individuals due to immunity during the latent period k, the rate at which exposed individuals become infected after the latent period α and the natural death rate. In turn, a proportion α of the exposed individuals become infected.

On the other hand, the rate of change of infected individuals decreases in a proportion equal to r + µ + d 1 , where r is the recovery rate of exposed individuals due to immunity and d 1 represents the death rate of infected individuals due to COVID-19. Also, an increase in the speed of the amount of infected individuals is considered in view that a proportion α of exposed individuals become infected. Meanwhile, the speed at which people are quarantined is increase by a fraction q 1 of exposed individuals, and decreases at a a rate equal to (q + µ + d 2 )Q, where q is the rate of recovery from quarantine to recovered and d 2 is the death rate of quarantined individuals due to COVID-19. Finally, the rate of change of the recovered individuals increases with the amounts of exposed, infected and quarantined who recover, and decreases with the portion of recovered members of the population who die due to natural causes. For convenience, Table 1 provides the set of model parameters employed and their physical relevance.

Using the assumptions in the paragraphs above, the mathematical model proposed to describe the dynamics of propagation of COVID-16 is given by the nonlinear system of coupled ordinary differential equations

Natural death rate of individuals k

Recovery rate of exposed individuals due to immunity during latent period r

Recovery rate of exposed individuals due to immunity α Rate at which exposed individuals become infected after latent period d 1 Death rate of infected individuals due to COVID-19 d 2

Death rate of quarantined individuals due to COVID-19 β 1

Contact rate of susceptible individuals with infected individuals β 2 Contact rate of susceptible individuals with exposed individuals q 1 Quarantined rate of exposed individuals q 2 Quarantined rate of infected individuals q Rate of recovery from quarantine to recovered 

and R 0 be nonnegative real numbers which physically represent the initial sizes of the susceptible, exposed, infected, quarantined and recovered subpopulations, respectively. More precisely, the following conditions are satisfied:

It is easy to see that, in the absence of COVID-19, the system

Notice that the reproductive number R 0 is the spectral radius of FV −1 , where

(2.5)

More precisely, notice that

.

(2.6)

Using the possible values of the reproductive number, we will devote our next efforts to establish conditions that guarantee the local asymptotic stability of the equilibrium points obtained above. Beforehand, notice that the Jacobian matrix associated to (2.1) is given by

For the remainder, the symbol Ω will represent the set of all (S , E, I, Q, R) ∈ R 5 such that S , E, I, Q, R ∈ R + . Theorem 1. The equilibrium point C 1 is locally asymptotically stable if and only if R 0 < 1.

Proof. Using the Jacobian matrix (2.7) and an easy substitution, we can readily obtain that

(2.8)

Notice the that the characteristic polynomial is given by

where a 1 = λ/µ, a 2 = κ+α+µ+q 1 and a 3 = r+µ+d 1 . The stability when R 0 < 1 follows now from the Routh-Hurwitz criterion for second-order polynomials, noticing that a 3 + a 2 − β 2 a 1 > 0 and (κ + α + µ + q 1 ) + (r + µ + d 1 ) − β 2 λ µ > 0 and, thus, all the eigenvalues are negative numbers.

Proceeding as in the proof of Theorem 1, it is easy to check that the characteristic polynomial is given by

where

(2.14)

(2. 16) Assuming that R 0 > 1 then the inequalities c 1 > 0, c 3 > 0 and c 1 c 2 > c 3 are all satisfied. It follows that all the eigenvalues of the problem (2.1) are negative. The conclusion of this theorem readily follows now using the Routh-Hurwitz criterion for third-order polynomials.

In the following theorems, we tackle the global asymptotic stability of the equilibrium points.

The equilibrium point C 1 is globally asymptotically stable if and only if R 0 < 1.

Proof. Consider the Volterra-type Lyapunov functional V : Ω → R, given by

(2.17)

A straightforward derivation shows that

(2.18)

Notice that the right-hand side of (2.18) is non-positive only if R 0 < 1, and it is equal to zero if S = S 0 * and E = I = Q = R = 0. This means that the only trajectory of the system on which V = 0 is C 1 . LaSalle's invariance principle assures then that C 1 is globally asymptotically stable in Ω.

If R 0 > 1 then C 2 is globally asymptotically stable, and it is unstable otherwise.

Proof. Consider now the Volterra-type Lyapunov functional U : Ω → R defined by

(2. 19) Notice now that

(2.20)

The right-hand side of (2.20) is obviously non-positive if R 0 > 1, and it is equal to zero only if S = S 1 * , E = E 1 * ,

This implies that the only set in which U = 0 is C 2 . LaSalle's invariance principle implies now that C 2 is globally asymptotically stable in Ω.

Before closing this section, we examine the sensitivity of the reproduction number with respect to each of the parameters involved. To that end, the following identities can be easily verified:

.

(2.22)

Also,

.

(2.25)

Observe that the numbers A β 1 , A β 2 and A λ are positive. Meanwhile, the remaining numbers are negative if R 0 < 1. We conclude that the sensitive parameters of the reproduction number are β 1 , β 2 and λ.

For each N ∈ N, define the sets I N = {1, . . . , N} and I N = I N ∪ {0}. The present section will be devoted to provide and analyze a nonstandard finite-difference discretization of the model (2.1). To that end, we let N ∈ N and consider a temporal period T > 0. Fix a uniform partition of the temporal interval [0, T ] consisting of N subintervals, and let τ = T/N. Define t n = nτ, for each n ∈ I N . For each n ∈ I N , we will employ the notations S n , E n , I n , Q n and R n to represent numerical approximations to the values of the functions S , E, I, Q and R, respectively, at the time t n .

Let U be any of the functions S , E, I, Q or R. We define

Using this nomenclature, the non-local finite-difference discretization of the mathematical model (2.1) is given by

δ t R n = κE n + rI n + qQ n − µR n+1 , ∀n ∈ I N−1 .

As expected, we will employ the discrete initial data (S 0 , E 0 , I 0 , Q 0 , R 0 ), where

Our first result establishes that the numerical method (3.2) is a positivity-preserving technique. Moreover, the explicit character of the scheme is exhibited in the proof of that result. Theorem 6. If the initial data satisfy (S 0 , E 0 , I 0 , Q 0 , R 0 ) ∈ Ω then the system (3.2) is uniquely solvable, and the solutions belong in Ω.

Proof. Beforehand, notice that the system (3.2) can be rewritten equivalently as

, ∀n ∈ I N−1 ,

, ∀n ∈ I N−1 ,

, ∀n ∈ I N−1 , R n+1 = R n + τ(κE n + rI n + qQ n ) 1 + τµ , ∀n ∈ I N−1 .

Observe that (S n+1 , E n+1 , I n+1 , Q n+1 , R n+1 ) is uniquely determined and belongs to Ω in the case that (S n , E n , I n , Q n , R n ) belongs to Ω. The result follows now by induction using the hypothesis on the initial data.

Next, we prove that the stability character of the equilibrium point C 1 is preserved by the method (3.2). The following technical result will be of utmost importance to that end. Lemma 7 (Brauer et al. [31] ). Let P 1 , P 2 ∈ R, and suppose that λ 1 , λ 2 ∈ C are the roots of the λ 2 − P 1 λ + P 2 . Then |λ i | < 1 for i = 1, 2 if and only if the following conditions are satisfied:

Proof. The proof of this result is straightforward.

(a) If R 0 < 1 then C 1 is locally asymptotically stable. (b) If R 0 > 1 then C 2 is locally asymptotically stable.

Proof. We will prove only (a) in view that the proof of (b) is similar though more tedious. Consider the right-hand sides of the equations in (3.2) as functions of S n , E n , I n , Q n and R n . Calculate then the Jacobian matrix and evaluate it when (S n , E n , I n , Q n , R n ) = C 1 . In such way, we obtain that

(3.5)

In turn, the characteristic polynomial of J is given by

where λ 1 = λ 2 = (1 + τµ) −1 , λ 3 = [1 + τ(q + µ + d 2 )] −1 . Moreover,

, (3.7)

. (3.8) Clearly, the absolute value of the eigenvalues λ 1 , λ 2 and λ 3 is less than one. To check that the the roots of the polynomial λ 2 − P 1 λ + P 2 also have absolute values less than one, it suffices to show that the conditions (i)-(iii) of Lemma 7 hold. However, the fact that P 1 > 0 assures that (i) is satisfied when (ii) of that lemma is true. To demonstrate (ii) and (iii), various simple algebraic calculations need to be performed, and the condition R 0 < 1 needs to be employed in both cases. As a consequence, all the eigenvalues of J have absolute values less than 1. This implies that C 1 is locally asymptotically stable, as desired.

Before establishing the main numerical properties of the scheme (3.2), we would like to plot the largest eigenvalues of the Jacobian matrix associated to C 2 . To that end, we will employ the estimated and fitted values of the parameters recorded in Table 2 . The graph of the spectral radius of the Jacobian matrix is plotted in Figure 2 as a function of the temporal step size. Notice that the spectral radius is always less than one, as predicted by Theorem 8. Table 2 were used, and the spectral radius was obtained for various temporal step-sizes τ in [0, 1000].

Proof. Beforehand, let us define the differential operators

On the other hand, define also the discrete operators

(3.10)

Here, we used the discrete operator

where U can be any of the functions S , E, I, Q or R for the remainder of this proof. Notice that, under the regularity assumptions and using the classical argument based on Taylor's theorem, it follows that there exist constants C 1 U and C 2 U which are independent of τ, such that

Using the smoothness of E and I, we can fix respective bounds C 0 E and C 0 I of these functions on [0, T ]. Using the triangle ineqlauity and the bounds above, we observe then that there exists a constant C U which is independent of τ, such that the following are satisfied for U being any of S , E, I, Q or R:

(3.14)

To be more precise, those constants are

We conclude that (3.2) yields first-order approximations to the solutions of (2.1), as desired.

The following is a discrete version of Gronwall's inequality.

Lemma 10 (Pen-Yu [33] ). Let (ω n ) N n=0 and (ρ n ) N n=0 be finite sequences of nonnegative mesh functions, and suppose that there exists C ≥ 0 such that

(3.20)

Then ω n ≤ ρ n e Cnτ for each n ∈ I N .

Theorem 11. There exists a τ 0 > 0 such that the finite-difference method (3. 2) is stable for values of τ < τ 0 .

Proof. Consider two sets of initial conditions (S 0 , E 0 , I 0 , Q 0 , R 0 ) and ( S 0 , E 0 , I 0 , Q 0 , R 0 ). The corresponding numerical solutions obtained through (3.2) will be denoted by (S , E, I, Q, R) and ( S , E, I, Q, R), respectively. Obviously, the first solution satisfies (3.2) while the second satisfies the discrete problem

δ t R n = κ E n + r I n + q Q n − µ R n+1 , ∀n ∈ I N−1 .

For convenience, we will define ε n U = U n − U n and ε 0 U = U 0 − U 0 , for each ∈ I N and U = S , E, I, Q, R. Calculating the difference between the respective equations of (3.2) and (3.21) , it is easy to see that the following problem holds:

(3.22)

Using the mean value theorem, it is easy to show that there exist constants C n S ,1 , C n S ,2 , C n S ,3 , C n E,1 , C n E,2 , C n E,3 ∈ R, for which the following identities are satisfied:

(β 1 I n + β 2 E n )S n+1 − (β 1 I n + β 2 E n ) S n+1 = C n S ,1 ε n I + C n S ,2 ε n E + C n S ,3 ε n+1 S , ∀n ∈ I N−1 , (3.23) (β 1 I n + β 2 E n )S n − (β 1 I n + β 2 E n ) S n = D n S ,1 ε n I + D n S ,2 ε n E + D n S ,3 ε n S , ∀n ∈ I N−1 , (3.24)

Let us consider firstly the first equation of the system (3.22) , and multiply both sides by τ. Take then the absolute value on both sides and use the triangle inequality on both. In that way, we readily obtain |ε n+1 S | − |ε n S | ≤ C n S τ(|ε n I | + |ε n E | + |ε n+1 S |) + µτ|ε n+1 S |, ∀n ∈ I N−1 .

(3.25)

Here, we let C n S = max{C n S ,1 , C n S ,2 , C n S ,3 }. Let now k ∈ I N−1 , and take the sum on both sides of the last inequality over all n between 0 and k. Using the formula for telescoping sums and letting C S = µ + max{C n S : n ∈ I N }, we reach

As a consequence, we reach that

In similar fashion, the remaining identities of (3.22) lead to the following inequalities

Let τ 0 = min{C −1 U : U = S , E, I, Q, R}, and let 0 < τ < τ 0 . Obviously, the numbers 1 − C U τ are positive for all U = S , E, I, Q, R, and we let C * be the minimum of all of them. Moreover, let C * = max{C U : U = S , E, I, Q, R}, and agree that ε n = |ε n S | + |ε n E | + |ε n I | + |ε n Q | + |ε n R |, for each n ∈ I N . Using then the inequalities above, it is easy to see that

ε n , ∀k ∈ I N−1 .

(3.32)

Using now Lemma 10, it follows that ε n ≤ Kε 0 , where K = C −1 * exp(C * T/C * ). As a consequence, the conclusion of the theorem readily follows now. Proof. The proof is similar to that of Theorem 11.

Before closing this section, we employed the finite-difference scheme (3.2) to confirm the validity of the analytical results. Figure 3 shows the solution of the system (2.1) using the method proposed in this work. To that end, we employed the parameters of Table 1 , and the initial conditions are. Two cases were considered, namely, (a) a system free of COVID-19 and (b) a system with the presence of COVID-19. The graphs show that the solutions tend to the prescribed equilibrium points C 1 and C 2 , respectively. It is worth pointing out that these results are in full agreement with both the analytical and the numerical theorems derived in this work. Figure 3 

In this work, we considered a generic SEIQR system to model the propagation of COVID-19. The model employed here consisted of a system of nonlinear ordinary differential equations of first order. Various constant parameters were used to describe the system and nonnegative initial conditions were employed. The method considered the presence of susceptible, exposed, infected, quarantined and recovered populations in light of the known dynamics of COVID-19. The mathematical model was obtained from epidemiological assumptions, and equilibrium points were obtained considering the presence and the absence of the epidemic. The reproductive number was obtained, and the stability of the equilibria was established in terms of the values of the reproductive number. Both local and global stability were analyzed. Also, a finite-difference method was proposed to investigate computationally the dynamics of COVID-19. The method is nonlocal, and we proved the positivity of the scheme, along with the discrete stability of the equilibrium points. In addition, the numerical properties of consistency, stability and convergence were rigorously established. Some numerical simulations proved the validity of the theorems on the stability of the equilibrium points. In this way, the present work contributes to the design of SEIQR models to describe the dynamics of infectious diseases like COVID-19 [34] .

The continuing 2019-nCoV epidemicthreatof novel coronaviruses to global health: the latest 2019 novel coronavirus outbreak in Wuhan, China

Handbook of COVID-19 prevention and treatment, The First Affiliated Hospital

Care for critically ill patients with COVID-19

COVID-19 and thrombotic or thromboembolic disease: implications for prevention, antithrombotic therapy, and follow-up

Keeping the lights on: Economic medicine for a medical shock

COVID-19 economic cost; impact on forcibly displaced people

Dynamics of a novel nonlinear stochastic SIS epidemic model with double epidemic hypothesis

Analysis of SIR epidemic model with information spreading of awareness

Traveling waves for a diffusive SEIR epidemic model with non-local reaction

Analysis and numerical simulations of a stochastic SEIQR epidemic system with quarantine-adjusted incidence and imperfect vaccination, Computational and mathematical methods in medicine

Numerical modeling and theoretical analysis of a nonlinear advection-reaction epidemic system

Analysis and nonstandard numerical design of a discrete three-dimensional hepatitis b epidemic model

A numerical efficient splitting method for the solution of two dimensional susceptible infected recovered epidemic model of whooping cough dynamics: Applications in bio-medical engineering

A deterministic model for the distribution of the stopping time in a stochastic equation and its numerical solution

Understanding epidemics from mathematical models: Details of the 2010 dengue epidemic in

An individual-based modeling framework for infectious disease spreading in clustered complex networks

Investigating spatiotemporal dynamics and synchrony of influenza epidemics in Australia: An agent-based modelling approach

Nonlinear dynamical analysis and control strategies of a network-based sis epidemic model with time delay

A stochastic SIRS epidemic model with nonlinear incidence rate

Some estimation problems in epidemic modeling

A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis

Positivity and boundedness preserving numerical algorithm for the solution of fractional nonlinear epidemic model of HIV/AIDS transmission

Positivity preserving operator splitting nonstandard finite difference methods for seir reaction diffusion model

A novel time efficient structure-preserving splitting method for the solution of two-dimensional reaction-diffusion systems

A positive and bounded finite element approximation of the generalized Burgers-Huxley equation

Isolation, quarantine, social distancing and community containment: pivotal role for old-style public health measures in the novel coronavirus (2019-ncov) outbreak

A positivity-preserving nonstandard finite difference scheme for the damped wave equation

A new positivity-preserving nonstandard finite difference scheme for the DWE

Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations

Mathematical models in population biology and epidemiology

Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative

Numerical methods for incompressible viscous flow

Modelling the impact of COVID-19 in Australia to inform transmission reducing measures and health system preparedness

Acknowledgments. The corresponding author (J.E.M.-D.) wishes to acknowledge the financial support from the National Council for Science and Technology of Mexico (CONACYT) through grant A1-S-45928.