Forward invariant set preservation in discrete dynamical systems and numerical schemes for ODEs: application in biosciences

We present two results on the analysis of discrete dynamical systems and finite difference discretizations of continuous dynamical systems, which preserve their dynamics and essential properties. The first result provides a sufficient condition for forward invariance of a set under discrete dynamical systems of specific type, namely time-reversible ones. The condition involves only the boundary of the set. It is a discrete analog of the widely used tangent condition for continuous systems (viz. the vector field points either inwards or is tangent to the boundary of the set). The second result is nonstandard finite difference (NSFD) scheme for dynamical systems defined by systems of ordinary differential equations. The NSFD scheme preserves the hyperbolic equilibria of the continuous system as well as their stability. Further, the scheme is time reversible and, through the first result, inherits from the continuous model the forward invariance of the domain. We show that the scheme is of second order, thereby solving a pending problem on the construction of higher-order nonstandard schemes without spurious solutions. It is shown that the new scheme applies directly for mass action-based models of biological and chemical processes. The application of these results, including some numerical simulations for invariant sets, is exemplified on a general Susceptible-Infective-Recovered/Removed (SIR)-type epidemiological model, which may have arbitrary large number of infective or recovered/removed compartments.


Introduction
Dynamical systems are pervasive in the modelling of naturally occurring phenomena [25].
While there are many results of the theory of continuous dynamical systems that permit to model the features of real-life processes, capturing them with discrete dynamical systems can be challenging.One typical and fundamental example is the forward invariance of a closed set D contained in open set ⊂ R n with respect to the following system of autonomous ordinary differential equations: in which the independent variable is the time t ≥ 0 and the function f is smooth enough.
For convenience, we recall that the said forward invariant property means that S(t)x 0 ∈ D at all time t ∈ [0, ∞) and for any where S(t) denotes the solution operator and we assume that there exists a unique global solution S(t)x 0 for any initial condition x(0) = x 0 ∈ D.
A short account on the invariance problem is given in [27].The basic hypothesis for the invariance property ( 2) is a tangent condition that roughly states that at a point x on the boundary ∂D of D, i.e., x ∈ ∂D, the vector f (x) is either tangent to D or points into the interior, • D, of D. In 1942, Nagumo [21] formulated the tangent condition in the following form, using the distance d(x, D) from a point x to the set D: It took nearly three decades for the invariance problem to be revitalized, starting with the seminal works of Brezis [8], who used the definition (3), and that of Bony [7], who came up with the following subtler formulation of the tangent condition based on the inner product, x, y , in R n : ν(x), f (x) ≤ 0 for x ∈ ∂D, ( 4 ) where ν(x) is any outer normal to D at x (in the sense of Bony).Let us recall that a vector ν(x) = 0 is called outer normal vector to D at x ∈ ∂D if the open ball with center x + ν(x) and radius |ν(x)| has no intersection with • D. The straight line through x perpendicular to ν(x) is called a tangent.This general definition of normal vector and tangent does not require any smoothness of the boundary ∂D, which is rather convenient in applications.Let us note that the normal vector ν(x), similarly the associated tangent, need not exist and, if existing, need not be unique.However, ν(x) at the point x ∈ ∂D reduces to the classical outer normal if the boundary is smooth, so the tangent exists in the classical sense at x.
In the late sixties, there was a great deal of solving the invariance problem (and /or related issues of existence of solution of the initial value problem) under either formulation of the tangent condition (3) or (4), as seen from the contributions in [10,13,23], where the works of Bony and Brezis [7,8] were revisited, compared, and extended.Since then, the tangent condition has abundantly been used, particularly in mathematical biology, to show the forward invariant structure of the positive cone R n + with respect to ordinary differential equation models and even with respect to reaction diffusion models (see for instance [24,28]).However, to the authors' best knowledge, the use of the tangent condition to establish the forward invariant nature of sets under discrete dynamical systems has not been explored.
The purpose of this work is primarily to fill this gap.To this end we consider a discrete dynamical system where in the map, F ≡ F(h, x) ≡ F(h)(x) : R × R n → R n , the parameter h > 0 represents the time step size between two consecutive states x k and x k+1 .This map is assumed to be as smooth as needed in the two arguments.The specific class of discrete dynamical systems we will focus on is defined as follows [12]: 5) is called symmetric or (time) reversible if More generally, the discrete dynamical system is said to be reversible for all x whenever y = F(h, x) implies that x = F(-h, y), which means that We introduce a discrete analog of the tangent condition that in essence states the following: there exists h > 0 such that for all x ∈ ∂D and h ∈ (0, h], the vector F(-h, x)x does not point inside • D. Under this condition, we show that the set D is invariant under the map F(h, .)for every h in (0, h]. A natural question of interest is to link the above-stated discrete tangent condition to continuous dynamical systems.This brings us to the second purpose of this paper, which we achieve by considering the sequence (x k ) obtained recursively through Eq. ( 5) as approximations at the discrete times t k = hk, h > 0, k = 0, 1, 2, . . . of the solutions x(t) = S(t)x 0 of the continuous system (1).In this case, it is essential to assume that F is at least continuous on x and continuously differentiable on h, ( 7 ) and that it satisfies the consistency requirements, Our specific aim is that the numerical method ( 5) is reliable in the three important directions below.The numerical method must satisfy the discrete tangent condition whenever the tangent condition holds for the continuous system, thereby showing that the scheme preserves the forward invariant structure.Furthermore, the numerical method must be elementarily stable and second-order convergent.
Regarding the latter targeted property, we stress that the construction of higher-order nonstandard finite difference (NSFD) schemes that are dynamically consistent with the underlying features of the continuous differential equations, particularly for those with transient dynamics, is a pending problem.Several authors have attempted to address this problem.These include the nonstandard versions of the classical θ -method, with θ = 1/2, developed in [1][2][3]16], where the positivity of the discrete solution is achieved by also using Mickens' rule on a nontrivial denominator function for the discrete derivative [18,19].In the same vein, we mention the recent work [15], valid for a scalar differential equation, where the second-order accuracy and elementary stability are achieved by a modified nonstandard θ method, 0 ≤ θ ≤ 1, in which the nontrivial denominator function varies at each iteration.
The rest of the paper is organized as follows.In the next section, we state precisely the discrete tangent condition in two equivalent forms and establish the forward invariance property under time-reversible flows.We devote Sect. 3 to the construction of a NSFD scheme that is time-reversible, and show that, in essence, this property makes our scheme of second-order accuracy and elementarily stable.In the same section, it is further shown via the discrete tangent condition that the NSFD scheme preserves the forward invariance structure of the continuous dynamical system.The case of mass action models of biological and chemical processes is considered in detail in Sect.4, followed by Sect. 5 on a general SIR-type epidemiological model with arbitrary number of infective or recoverred/removed compartments.Concluding remarks with possible extensions of this work are given in Sect.6.

Domain preserving property of reversible schemes
The theory for continuous dynamical systems of the form (1) provides theorems proving that a set in R n is forward invariant through conditions only on the boundary of the set, namely the tangent condition as stated in (3) or (4), see for instance [27].In this section we derive discrete counterparts of the tangent condition and resulting theorems for the forward invariance of sets under reversible maps.
Let the closed set D be a subset of and be related to its interior, • D, and boundary, ∂D, as follows: Theorem 2 Let the discrete dynamical system (5) be reversible for all h ∈ (0, h] and let (7)-( 8) hold.Then the set D is invariant under F(h, •) for every h ∈ (0, h) if and only if, Proof To show that Condition (10) is sufficient for the invariance of D, let us consider ) is continuous on h, there exists ĥ ∈ (0, h] such that y = F( ĥ, x) ∈ ∂D.Equivalently, x = F(-ĥ, y), which contradicts the Condition (10).Therefore, (9) and the continuity of F on x, we obtain Conversely, let D be invariant under F(h, •) for every h ∈ (0, h], i.e., F(h, D) ⊆ D for every h ∈ (0, h], (11) and assume that Condition (10) is violated, that is, there exists h ∈ (0, h] and x ∈ F(-h, ∂D) ∩ • D.Then, using the fact that F(-h, •) is the inverse of F(h, •) we have From the continuity of F(-h, •), it follows that the set F(h, This is a contradiction to (11), which proves that Condition (10) is necessary for the stated invariance of D.
Condition (10) is not easy to verify.If the set D is convex, we replace below Condition (10) by an equivalent condition, requiring that the vector F(-h, x)x be outward directed or tangential to D at x ∈ ∂D.The latter condition is our proposed discrete analog of Bony [7] tangent condition for continuous dynamical systems, based on Bony's definition of outer normal vector ν(x), stated in (4).
Let D be convex.Then some comments on outer normal vectors are in order.First, let us recall that a supporting hyperplane to D at a point x ∈ ∂D is a hyperplane containing the point x, while the set D is contained in one of the closed halfspaces determined by the hyperplane.The so-called supporting hyperplane theorem [17] states that for D, a closed convex set with nonempty interior, there exists a supporting hyperplane at every x ∈ ∂D.It is easy to see that the vector ν(x), which is normal to the supporting hyperplane at x and points in the halfspace not containing D, is an outer normal vector to D at x in Bony's sense given in the introduction.Under the stated assumptions, it follows from the hyperplane separation theorem for convex sets [17] that one can associate a supporting hyperplane with every outer normal vector, a fact that is intuitively clear.To fix the notation, the supporting hyperplane at x ∈ ∂D with an outer normal vector ν(x) is the set which is the common boundary of the following two closed halfspaces: Then, the fact that or, equivalently, Theorem 3 Let a discrete dynamical system (5) be reversible for all h ∈ (0, h] and let D be a closed and convex subset of satisfying (9).Then the set D is invariant under F(h, •) for every h ∈ (0, h] if, and only if, for every h ∈ (0, h] and x ∈ ∂D, there exists y ∈ ∂D and an outer normal vector ν(y) such that Proof We show that the discrete tangent condition (14) implies Condition (10).Indeed, let h ∈ (0, h] and x ∈ ∂D.Using the convexity of D, the inequality (14) shows that Using ( 13), we obtain Then ( 15) and ( 16) Considering that h ∈ (0, h] and x ∈ ∂D are arbitrary, Condition (10) follows.
To show the converse implication ( 10) ⇒ ( 14), let h ∈ (0, h] and x ∈ ∂D.Condition (10) shows that From the hyperplane separation theorem of convex sets it follows that there exists y ∈ ∂D and an outer normal vector ν(y) such that F(-h, x) ∈ H + ν(y) , which implies (14).
The difficulty raised about the practical use of Condition (10) is equally present in its equivalent formulation (14).In applications, it is easier to use the latter condition with y = x, i.e., Eq. ( 17) below.This turns out to be our discrete tangent condition, which, as stated in the next theorem, is a sufficient condition under which the forward invariance property with respect to discrete dynamical systems is achieved.
Theorem 4 Let a discrete dynamical system (5) be reversible for all h ∈ (0, h] and let D be a closed and convex subset of satisfying (9).

Reversible nonstandard finite difference schemes
In the construction of nonstandard finite difference schemes, Mickens [18,19] made an important observation, namely, that the discrete models of differential equations have a larger parameter space than the corresponding differential equations.Adding to this Mickens' Rule 1, which says that the orders of the discrete derivatives must be exactly equal to the orders of the corresponding derivatives of the differential equations, it is convenient in the setting of Mickens' nonlocal approximation Rule 3, to view the right-hand side of the system (1) as a restriction, on the diagonal z = y, of a certain function of two variables ϕ ≡ ϕ(y, z) such that f (x) = ϕ(x, x) and hence the system takes the form We consider the numerical method, which is a NSFD scheme [4], on which only the nonlocal approximation of the right-hand side of ( 18) is used, the rule on the complex denominator function of the derivatives being excluded.The NSFD scheme ( 19) is implicit.The existence and uniqueness of solution of the respective system of algebraic equation is an issue which requires attention on its own.To make this NSFD scheme a generalized dynamical system on [25], we will assume here that Finding general conditions for (20) to hold is challenging, as is the case for general equations of the form G(x, y) = 0 solved by using for instance the implicit function theorem, sophisticated fixed-point iterations, etc.However, in particular settings, as the ones that arise in applications and are considered in the sequel, such conditions are relatively easy to formulate.
Theorem 5 Assume that (20) holds.Then the NSFD scheme (19) satisfies the following properties: a) It is time-reversible; b) It is elementarily stable in the sense that b1) Its fixed points are precisely the equilibrium points of the system (18); b2) It replicates the local stability of all hyperbolic equilibrium points of (18); c) It is a second-order scheme for the Equation (18).
Proof The Assumption (20) implies that the NSFD scheme (19) can be written in the explicit form (5), where the function F satisfies the equation a) If in the Equation ( 19) we interchange x k and x k+1 and replace h by -h, the equation remains the same.Therefore, (6) holds, which means that the NSFD scheme is timereversible.
b1) It is clear that x is a fixed point of ( 19) if, and only if, ϕ(x, x) = 0, i.e. x is an equilibrium point of (18).b2) Let us use the explicit form (5) of (19).From (21), the Jacobian matrix (in the x variable) of F satisfies where I denotes here and after the n × n identity matrix.Let x be an arbitrary fixed point of (19).At x = x, using F(h, x) = x, Equation ( 22) simplifies to or, equivalently, Let λ ≡ (λ) + ı (λ) be an eigenvalue of ∂f (x) ∂x with associated left eigenvector v. Then multiplying both sides of ( 23) on the left by v, we obtain It follows from ( 20) that the matrix ∂F(h,x) ∂x is not singular.Then, due to the obvious impossibility for 1 -h 2 λ and 1 + h 2 λ to be both zero, neither of them is.Hence Therefore, to every eigenvalue λ of ∂f (x) ∂x with associated left eigenvector v corresponds an eigenvalue μ = of the matrix ∂F(h,x) ∂x with the same left eigenvector v.After some technical manipulation, we obtain It follows from ( 24) that which implies that a hyperbolic equilibrium point x of ( 18) is asymptotically stable if and only if it is asymptotically stable as a fixed-point of (19).c) Let x ∈ and let u(t) = S(t)x denote the solution of (1) with u(0) = x.Denote Taylor expansion about h = 0 yields by the consistency conditions (8).Then, for the local truncation error, of the NSFD scheme ( 5) and ( 21), Taylor expansions of u about h 2 and of ϕ about (y, z) = (ξ , ξ ) yield Then performing a further Taylor expansion and using (25), we obtain E(h) = O(h 3 ) from which second-order accuracy follows by the standard theory for one-step numerical methods.
Remark 6 Since the scheme ( 19) is symmetric, the statement c) in Theorem 5 can be derived from the general theory of symmetric schemes [12].We note that properties in items a) and b) are of qualitative nature and do not follow from the standard numerical analysis theory, involving zero-stability, consistency, and therefore convergence.

Application to mass action-type models
Modelling of biological and chemical processes by applying the mass action principle for representing the interaction of the involved species typically results in a system of the form (18) in , a convex and compact subset of R n + , where the function ϕ is linear in both its arguments.Then ϕ can be represented as where P and Q are n × n matrix functions of x and y, respectively.The matrix A and the vector b are constant.In this particular setting, the numerical method ( 19) can be written in the form or, equivalently, Furthermore, and also expressing its reversibility, the equation ( 27) can be written in the form Considering that is compact, there exists h such that the matrices are both diagonally dominant matrices for all h ∈ (0, h] and x ∈ . (30) Then, these matrices are non-singular and ( 28) can be written explicitly as where Further, we have Theorem 3 is a useful tool for determining the invariance of under F(h, •).Using the specific form of the maps (32), we derive an explicit representation of F(-h, x)x.Indeed, using (33) and the expression of M(h, x) in (30), we have Therefore, and also, by Taylor expansion, Then, the relation ( 17) can be written equivalently as Assume that the continuous tangent condition (4) holds with strict inequality < 0 for all x ∈ ∂ and for the dynamical system (1) or (18) in which the function f (x) = ϕ(x, x) is given by (26).To be more precise, we assume that there exists a constant upper bound c < 0 in the said strict inequality, i.e., ≤ c < 0 for all x ∈ ∂ .Then, it follows from (37) that there exists h > 0 such that the relation (36) holds in the following uniform manner, which is the discrete tangent condition (17): On the contrary, if (4) holds in the limit case with the equality '= 0' , the discrete tangent condition (38) must be checked directly, because the previous argument is no longer true.Note that the scheme ( 27) is constructed in the general form (19).It satisfies Assumption (20) with the unique solution being given in the explicit form (31) for h ∈ (0, h].The invariance of follows from Theorem 3 combined with the discrete tangent condition (38), which, as shown above, is inherited from the continuous tangent condition in the strict case and is checked in the limit case.We have thus proved the following theorem: Theorem 7 Assume that the tangent condition (4) is satisfied for the continuous model (1), (18), and (26).Then the NSFD scheme (27) inherits the forward invariance on in the sense explained above.Furthermore, the scheme (27) satisfies the properties a)-c) in Theorem 5. Remark 8 In limit situations which are not conclusive (e.g., stability of the disease-free equilibrium when the basic reproduction number R 0 = 1 [2,14]), it is standard in mathematical analysis to deal with such cases differently.This is the essence of the requirement that the discrete tangent condition be checked directly in the limit case when ν(x), f (x) = 0.As a matter of fact, it often occurs in practical applications that the vector ν(x) is a left eigenvector of M(-h, x) with a positive eigenvalue λ.If this is the case, then the discrete tangent condition (38) follows directly from the tangent condition (4) of the continuous system.Indeed, Remark 9 The NSFD scheme ( 19) is an implicit method.When applied to mass action models, and due to the nonlocal approximation of nonlinear terms, it is interesting to observe that the method can be written in the form (28). Therefore, every step requires only solving a linear system-a computational effort comparable to the explicit methods.

Example: a general epidemiological model
In this section, we show how the theory and tools developed in Sects.2, 3, and 4 can be put together in deriving second-order scheme with the properties a)-b) in Theorem 5.
While the theoretical existence of h is important, deriving constructively the value of h is critical for any specific application.We show how the value of h can be computed for the considered model.Although the numerical approach proposed in this paper is exemplified on the model considered in this section, its realm of application is wider.

The model
We consider a general compartmental model of a contagious disease in a single population.Let the model have n compartments with their sizes given by the coordinates of a vector x = (x 1 , x 2 , . . ., x n ) T ∈ R n , where we fix x 1 = S-the compartment of susceptible individuals.
We assume that the infection rate is given via mass-action term of the form If the jth compartment contains infective individuals, then β j > 0. If the individuals in the jth compartment are not infective, then β j = 0. Examples of infective compartments are: symptomatic infective, asymptomatic infective, recovered, but still infective, partially isolated infective, etc. Examples of not infective, other than the susceptible, are: exposed (latent period), recovered with temporary immunity, recovered with permanent immunity, fully isolated (receiving treatment or not), dead, vaccinated (one or more types of vaccines), and others.The newly infected individuals are distributed in the compartments x 2 , . . ., x n .Let g ij Sx j be the growth rate of the ith compartment due to the susceptible infected by the individuals in the jth compartment.Clearly, we have g ij ≥ 0, i, j = 2, . . ., n, and Denote g 1j = -β j , j = 2, . . ., n, and g i1 = 0, i = 1, 2, . . ., n.Then using the matrix G = (g ij ), the mass action is represented by the vector x 1 Gx.The remaining transfer rates between compartments are assumed to be linear and represented by a matrix C, which is a Metzler matrix, that is, the nondiagonal entries of C are nonnegative.Further, a constant recruitment rate is given by a vector b = (b 1 , b 2 , . . ., b n ) T ≥ 0, where b i is the recruitment rate in the ith compartment.The recruitment in the compartment of susceptible is always positive, that is b 1 > 0, but other compartments may also have positive recruitment rates.
Then, the compartmental epidemiological model is represented as a system of differential equations where e 1 is the first vector in the canonical basis {e 1 , e 2 , . . ., e n } in R n .There are a large number of models in mathematical epidemiology that can be represented in the generic form (40).These include the SEIR-types of models, [14], as well as models with more strains, [9,26], models with multiple types of infective and treated individuals (SEITR model), [22], models including vaccination (SVEIR model), [11].
The system (40) is of the form (18) with which in turn is of the form (26), where and vector b as given.Assuming that the compartments take jointly into account the whole population, by adding all equations we obtain that the total population N = x 1 + x 2 + • • •+ x n satisfies the conservation law, and μ is the disease-independent removal/death rate.It is easy to see that this property implies that ν = (1, 1, . . ., 1) T is a left eigenvector of C and we have ν T C = -μν T .Since C is a Metzler matrix, this shows that all diagonal entries of C are negative.More precisely, c ii ≤ -μ < 0. Let us mention that due to the property (39), we have ν T G = (0, 0, . . ., 0) T .Hence, the mass action kinetics do not affect the demography.
It is easy to see, using the tangent condition ( 4), that the model (40) defines a forward dynamical system on Indeed, is bounded by the coordinate planes and the plane The respective outer normal vectors are -e 1 , -e 2 , . . ., -e n , ν.
Using the Metzler property of C, we have Further, on the plane (42), we have which completes the proof of the claim that is forward invariant for (40).

Discretization and its properties
We consider the numerical method (27) for the model (40) on .The matrices M(h, x) and M(-h, x), see (30), are both strictly diagonally dominant if and only if Therefore, if then both matrices M(h, x) and M(-h, x) are strictly diagonally dominant for all x ∈ and h ∈ (0, h].Hence the scheme can be written in the explicit form (31).
Our aim is to apply Theorem 7. Here, we will derive the invariance of under (31) directly since one of the tangent conditions is satisfied as an equality, see (43).In this way we also derive a computable value of h.
We deal first with the plane (42).Considering that ν T P = ν T Q = (0, . . ., 0) T and that Therefore, if h satisfies (44), ν is a left eigenvalue of M(-h, x) with a positive eigenvalue 1 -h 2 μ.Then the tangent condition (38) holds on the entire plane (42) due to Remark 8.
Next, we consider the boundary of on the coordinate plane x 1 = S = 0.It follows from (35) that the vector y = x -F(-h, x) is the solution y = (y 1 , . . ., y n ) of the linear system or, equivalently, For the tangent condition to hold, we need to show that Using that x ∈ and that y is a solution of the linear system (45), we obtain Assume that (44) holds.Then, (46) yields From the first equation in the system (45), we have Taking into account that the coefficient of y 1 is positive, y 1 is nonnegative provided Then we have established so far that the tangent condition (38) holds on the entire plane (42) and on the boundary of on the coordinate plane x 1 = 0 for all h ∈ (0, h].
The fact that F i (h, x) ≥ 0 for i = 2, . . ., n, x ∈ , h ∈ (0, h] can be established in a more direct way.Let the (n -1) × (n -1) matrix M(h, x) be obtained from M(h, x) given in (30) by deleting the first row and the first column.It is easy to see that M(h, x) = I -h 2 (x 1 G + C), where G and C are obtained from matrices G and C by deleting their first row and the first column, respectively.We note that the non-diagonal entries M(h, x) are non-positive.Further, M(h, x) inherits from M(x, h) its diagonal dominance.Hence, M(h, x) is an Mmatrix, which implies that ( for x ∈ and F 1 (h, x) ≥ 0. Let us recall that, while we established the tangent property for the whole plane (42), we established the tangent property only for that part of the coordinate plate x 1 = 0, which is boundary of .Hence, this result by itself does not imply F 1 (h, x) ≥ 0. We establish the nonnegativity of F(h, x) as follows.The tangent condition on the boundary of on x 1 = 0, implies that as h increases from 0 to h, F(h, x) may not leave through this part of the boundary.Let x 1 > 0. Then F 1 (h, x) > 0 for h small enough.It follows from (48) that, as h increases, while F 1 (h, x) > 0, we have F(h, x) ≥ 0. Hence, F(h, x) may leave only through the boundary on the plane x 1 = 0. Due to the tangent condition, this may not happen for h ∈ (0, h].Therefore, F(h, x) ∈ for all x ∈ and h ∈ (0, h]. Thus, from Theorem 7, we obtain that the NSFD scheme for the model ( 40) is reversible, it defines a discrete dynamical system on the domain of the continuous dynamical system (40), it is elementarily stable and second-order accurate.
Remark 10 The model (40) assumes only one compartment of susceptible individuals.In principle, the same method of proof works for a vector of susceptible compartments (for instance, as in [20]).We present this simple version of the model to avoid technicalities, which may obscure the main ideas demonstrated on this example.Further, we note that if there there is a disease induced mortality, (41) is satisfied as an inequality.In order to have (41) as equation, we can transform the model by introducing a compartment or compartments for deaths due to the disease with a removal rate equal to the natural death rate.This transformation increases the dimensions of the domain of the model, but allows for the case of disease induced mortality to be accommodated within the framework of ( 27) applied to the model (40).
The proof of the properties of NSFD method ( 27) applied to the model (40) uses the specific demographic equation (41).Let us note that demographic dynamics are not a restriction to applying the NSFD method (27).However, different demographic equations would require different approaches to proving the tangent condition for the domain of the model.
Essential for the application of the NSFD method (27) is the fact that the force of infection is represented through the mass action principle and is given in the form n j=2 β j x j .We may note that a model with standard incidence force of infection in the form n j=2 β j x j N , is easily converted to mass action model for z = x N , [9,14].Hence, the NSFD method ( 27) can be applied.However, the transformation changes the model and the domain, so that the proof of the tangent condition need to be adapted accordingly.Clearly, if the model cannot be written in a mass action form, the NSFD method (19) can be applied under the Assumption (20), but the method ( 27) is not applicable.
Remark 11 It follows from Theorem 7 that the NSFD scheme preserves all hyperbolic equilibria of the model as well as their stability properties, without any prior knowledge of the values of these equilibria or even their existence.In this way, the discetizations by (27) provide means of discovery of asymptotic dynamics of the model.This is solely done by using the nonlocal discretizations of nonlinear terms; a way to bring in our NSFD scheme Mikens' rule on nontrivial denominator function is discussed in [5].A natural question is if this result cannot be extended to invariant sets other than equilibria.There are nonstandard numerical methods preserving, for specific models, specific invariant sets such as first integrals considered for instance by the authors [4].However, constructing schemes and developing the associated theory, for preserving invariant sets in general is an open problem.Nevertheless, it is interesting that the method (27) applied to (40) preserves two important invariant sets on the boundary of the domain, as discussed in the next subsection.

Preserving invariant sets
The model (40) has two forward invariant sets on the boundary of the domain : the disease-free manifold and the constant population manifold.We show that both are preserved by the numerical method (27).
The disease free manifold is given by where I is the set of the indexes of the compartments promoting the disease-infective or transferring into infective.If the disease-free manifold exists, then the model (40) has the following properties if e j , x = 0, j ∈ I, then Gx = 0 and (e j , Cx = 0, j ∈ I. (51) Let x ∈ DFM.Then, Using (50)-( 51), it is easy to see that e j , I + h 2 C x + hb = 0, j ∈ I, and, consequently, by mathematical induction, Since h ≤ h, we have It follows from ( 52) and ( 53) that e j , F(h, x) = 0, j ∈ I. Considering that it has been proved already that F(h, x) ∈ , we obtain that The constant population manifold is given by Considering the specific form of the matrix M(h, x) for (40), multiplying (28) on the left by ν, we obtain or, equivalently, ) can be written as which is a second-order finite difference method for the Equation (41).

Some illustrative numerical simulations
In Sect.5.2 and Sect.5.3, we established that the NSFD scheme ( 28) applied to the model (40) preserves some properties of (40) considered as a dynamical system.These include: the domain , the hyperbolic equilibria with their stability, the disease-free manifold and the constant population manifold.These are properties of qualitative nature and cannot be derived through the techniques of the standard numerical analysis, involving consistence, stability and therefore convergence.To illustrate this point, we consider examples where well-known methods do not preserve the mentioned properties for a very simple version of (40), the SIR model: For the simulations, we use (i) midpoint method, which is a second-order method, and (ii) the fourth-order Runge-Kutta method as provided by the function ode45 in Matlab.Figures 1 and 2 provide phase diagrams of some numerical solutions in the S × I plane.The boundary of the domain is marked by a solid magenta line.The red star denotes the disease-free equilibrium, which is a stable equilibrium of the model (56) for the stated values of the parameters.Both simulations illustrate that numerical solutions produced by the stated methods may leave the domain of the model.Further, the trajectories by the Midpoint method presented in Fig. 1 converge (correctly) to the stable equilibrium, but some of them approach it from the exterior of the domain.The numerical solution presented in Fig. 2 eventually oscillates with nondecreasing amplitude as it approaches the stable equilibrium.Note that the Matlab function ode45 automatically selects optimal step size.Hence, h is not specified.
For comparison, we present on Fig. 3 the numerical solutions by the method ( 27) with h = 0.1 using the same data as in Fig. 1 and in Fig. 2. For the model (56), we have h = h = 2 γ μ +ν+δ .Hence, for both sets of values of parameters we have h = 0.1 < h.Therefore, it follows from the discussion in Sect.5.2 that the method (27) satisfies the properties stated in Theorem 7. One can also observe on the graphs in Fig. 3 that both the domain  and the stability of the equilibrium of the model (56) are preserved.Further, we may note that the values of the parameters used in Fig. 1 are associated with a disease, which is low contagious and with a fast removal rate.Hence, all solutions very quickly approach the disease-free manifold, which is the primary driver of the long-term dynamics.This property is correctly represented by the trajectories on Fig. 3(a), but not by the trajectories on Fig. 1, including those within the domain.

Conclusion
The literature on forward invariant sets in R n with respect to continuous dynamical systems defined by ordinary differential equations is rich (see, for instance, the short account in [27]).However, the invariance problem is not sufficiently investigated for discrete dynamical systems.Our findings in this work are threefold.Firstly, we formulated a discrete analog of the tangent condition by Bony [7], under which we established the forward invariance of sets with respect to time-reversible discrete dynamical systems.Secondly, we constructed a new NSFD scheme to approximate the solutions of continuous dynamical systems.Apart from its elementary stability, this new scheme is innovative and reliable in that: (a) it is time reversible, (b) it inherits in a discrete manner the tangent condition and the associated forward invariant property from the continuous system, and (c) it has second-order accuracy, thereby addressing the challenge of designing high-order numerical schemes that cannot exhibit spurious/ghost solutions or other elementary instability that do not correspond to the feature of the continuous equations.Thirdly, we showed that our construction and findings apply directly for mass action-type models in biological and chemical processes, specifically the preservation by the NSFD scheme of forward invariance of the biologically feasible domain and other dynamics of the continuous system.
Our future research includes the following: (1) Formulating a discrete analog of the tangent condition (3) by Nagumo [21] in the case when the latter is not equivalent to Bony's condition (4).Conditions for their equivalence are stated in [27].(2) Investigate the preservation by our time-reversible NSFD scheme of the global asymptotic stability of the disease-free equilibrium of continuous epidemic models when the basic reproduction number is less than 1. (See [2] for an approach, which avoids the use of the Lyapunov function for discrete dynamical systems.)(3) Extending this study to dynamical systems with non-hyperbolic equilibrium points, an approach considered in [6] as far as the stability for one-dimensional models is concerned.(4) Though a standard incidence-based epidemic model can be transformed into a mass action-based one (with different dynamics) to which the NSFD scheme proposed in this paper can be applied, as mentioned in Remark 10, investigating the reverse process in order to recover the preserved dynamics of the initial model is an issue of interest.(5) Other possible generalizations are: (i) NSFD scheme (27) applied to extended model (40) with vector of susceptible compartments and/or disease-induced death rate, (ii) models with demographics different from (40), and (iii) models where the general method (19) can be applied but not the method (27).

Figure 3
Figure 3 Phase diagrams of the numerical solutions by the NSFD method (27) of (56) with h = 0.1 and parameters as in (a) Fig. 1 and as in (b) Fig. 2 is open.Therefore, there exists ε > 0 such that the open ball V ε y with radius ε and centered at y satisfies V ε (y) ⊆ F(h, Since y is on the boundary of D, the open ball V ε (y) contains points which are not in D.