A non-standard discretized SIS model of epidemics

1.   Introduction

In most cases mathematical modeling of epidemic dynamics is based on formulating continuous systems. Discrete systems are less often applied. However, using computer software for simulations is becoming increasingly common. Furthermore, the epidemiological data reflect values of particular magnitudes only for determined moments. Because of these reasons, the discrete approach is gaining attention. The examples of discrete models describing epidemic phenomena can be found e.g., in ^[1,2].

Generally discrete models are constructed on the grounds of their continuous counterparts. One of the easiest method of discretization is the explicit Euler method (EEM). However, applying this method can lead to dynamical inconsistency of discrete models with their continuous analogues. In particular, there can arise biologically unreasonable bifurcating behavior. Moreover, negativity of system's solutions can occur. To avoid these nuisances, one can apply the non-standard discretization method (NSDM), which is presented and discussed in ^[3]. Discrete models exploiting this method are dynamically consistent with their continuous counterparts ^[4].

Even in systems discretized using the NSDM bifurcations can occur. Here we will focus on the possibility of appearing Neimark-Sacker bifurcation (NSB) in the analyzed system. We say that in a system of difference equations there can occur NSB if the Jacobian matrix for a given stationary state has two coupled complex eigenvalues with non-zero imaginary parts and modulus equal to $1$ ^[5]. This bifurcation is a discrete counterpart for Hopf bifurcation in continuous systems. However, the Hopf bifurcation is a term used also for a discrete case ^[6].

In NSB we identify a bifurcation parameter, which critical value gives desirable eigenvalues. There are two types of NSB: supercritical and subcritical one. In the supercritical NSB for bifurcation parameter's values below critical one a given stationary state is a sink. After exceeding this critical value this state loses stability and there arises a stable isolated closed invariant curve which surrounds the stationary state. In the subcritical bifurcation before exceeding critical value there is an analogous unstable curve surrounding a stable stationary state. After exceeding this value, the curve disappears and the stationary state loses stability ^[5].

Here we will consider a system belonging to the class of SIS (susceptible-infected-susceptible) models. In this kind of models it is assumed that there is no long-lasting immunity and immediately after infection individuals become susceptible again. The exemplary basic SIS models were presented in ^[7,8].

Authors in ^[9] analyze a discrete SIS model with the use of the EEM. The authors use a disease incidence function relied on the mass-action law and assume that growth of a population is constant. Dynamically consistent SIS models are presented in ^[10]. Here the standard incidence function was chosen. The discrete SIS model from ^[11] with the NSDM is analyzed in the context of the basic reproduction number. A similar approach is used in ^[12], where authors use a variation of the backward Euler discretization. Currently the main contribution of analysis of SIS models is related to epidemic spread in complex networks. In ^[13] authors investigate a continuous model for such a network, assuming that infection rates of multiple edges interfere with one another. A continuous model is analyzed also in ^[14], where an infective medium is considered. A discrete approach with the EEM is presented in ^[15]. However, in this three papers immigration, birth and death processes are neglected. Dynamics of the current COVID-19 pandemic is also an issue investigated from the mathematical point of view. In ^[16] both discrete and continuous models are introduced. Here the author includes also an additional group of recovered people, proposing SIR type (susceptible-infected-recovered) model. An incubation period is considered and an incidence function for a continuous case is built with the use of an integral.

In our previous papers (see ^[17,18]) we analyzed continuous SIS models differing in a type of inflow into subpopulations. Basing on the system from ^[17], we constructed the discrete model with the use of the EEM. This model was discussed in ^[19]. Here we will show a modification of this discrete model with the use of the NSDM.

This paper is organized as follows. In the next section we quote the continuous model from ^[17] and the discrete model with the use of the EEM from ^[19]. Then we introduce the non-standard discretized model and show the basic properties of this model. Later we focus on the stability analysis of stationary states of the system. Finally, we discuss the possibility of NSB appearance in the system. In the last part we give conclusions. Our theoretical results are complemented with numerical simulations with the use of the Matlab software.

In the paper we will assume that $\mathbb{N}$ means the set of all natural numbers including zero and $\mathbb{N}_+ : = \mathbb{N} \setminus \{ 0 \}$.

2.   A continuous model and a discretized model with the explicit Euler method

Let us first remind the continuous model analyzed in ^[17]. This model describes the epidemic spread in a homogeneous population in which we consider two exclusive groups: healthy (susceptible) and infected people. The sizes of these groups at time $t$ are denoted by $S(t)$ and $I(t)$, accordingly. Naturally, these sizes can be interpreted as densities. The model has the form:

where $C$ reflects an inflow into the population and coefficients: $\beta$, $\gamma$, $\mu$ and $\alpha$ correspond to illness transmission, recovery, natural death and disease-related death, respectively. The parameters appearing in system (2.1) are constant and positive. Because of the meaning of the parameters, it is reasonable to assume $C \gg \mu$.

We conducted scaling of system (2.1) in order to reduce the number of parameters. Eliminating $\beta$ and $\gamma$, we obtain

where $x = \beta S$, $y = \beta I$, $k = 1 + \alpha+ \mu$ and the coefficients in system (2.2) are scaled accordingly.

Let us remind the basic properties of system (2.2). The form of the right-hand side of system (2.2) guarantees existence, uniqueness and positiveness of solutions for positive initial conditions. The basic reproduction number $\mathcal{R}_0$ for system (2.2) is equal to

We will refer to this value also in the case of discrete systems.

System (2.2) has two stationary states:

1) disease-free: $E_d = (x_d, y_d) = \left (\tfrac{C}{\mu}, 0 \right)$, always existing,

2) endemic: $E_e = (x_e, y_e) = \left(k, \tfrac{C - \mu k}{k-1} \right)$ which exists for $C > \mu k$, that is for $\mathcal{R}_0 > 1$.

The state $E_d$ is locally stable for $\frac{C}{\mu} < k$. If $\frac{C}{\mu} > k$ (i.e., when $E_e$ exists), then $E_d$ is a saddle point. The state $E_e$ is always locally stable under its existence. Furthermore, the Poincaré Theorem implies global stability – either $E_d$ is globally stable whenever $E_e$ does not exist, or $E_e$ is globally stable whenever exists.

Now let us remind the discrete version of system (2.2) with the use of the EEM. This model was introduced and analyzed in ^[19]. Let $n \in \mathbb{N}$ be an $n$-th node in a discrete timescale. After discretization we obtained the system

where $h$ is a step size of discretization.

Now we quote the properties of system (2.4). Basing on the approach presented in ^[20] dedicated for discrete systems, we computed $\mathcal{R}_0$ for system (2.4) and obtained again system (2.3). Let us define a variable $w_n$ as

meaning the size of the whole population for the $n$-th point in the discrete timescale. Adding by sides the equations from system (2.4), we get

Solving inequality (Eq 2.5) for

yields

and hence we have

It means that the size of the whole population is bounded from above. Moreover, if inequality (Eq 2.6) holds and $w_0 \leq \frac{C}{\mu}$, then $w_n \leq \frac{C}{\mu}$ for every $n \in \mathbb{N}_+$.

Now let us remind the property related to the boundedness of the particular variables of system (2.4).

Proposition 1. If

and

then we have $x_n > 0$, $y_n > 0$ and $x_n+y_n \leq \tfrac{C}{\mu}$ for any $n \in \mathbb{N}_+$.

The forms of stationary states and conditions for their existence are the same in both systems (2.2) and (2.4). Let us present the conditions for stability of the stationary states in system (2.4).

Proposition 2. If inequality (Eq 2.9) holds, then $E_d$ is a sink for $\mathcal{R}_0 < 1$, non-hyperbolic for $\mathcal{R}_0 = 1$ and a saddle for $\mathcal{R}_0 > 1$. Moreover, if inequality (Eq 2.10) holds, then $E_d$ is globally stable for $\mathcal{R}_0 < 1$.

For the next proposition we define

and

Now we formulate the statement.

Proposition 3. Let us assume that $E_e$ exists. If $\delta \geq 0$, then $E_e$ is a sink for $h < h_a$, a saddle for $h_a < h < h_b$, a source for $h > h_b$ and non-hyperbolic for $h \in \{h_a, h_b \}$. If $\delta < 0$, then $E_e$ is a spiral sink for $h < h_c$, a spiral source $h > h_c$ and non-hyperbolic for $h = h_c$.

Comparing the conditions for the local stability of $E_e$ in system (2.2) to those in system (2.4), we state that the behavior of $E_e$ in the discrete system with the explicit Euler method is more complex than in the continuous system. Furthermore, in the discrete case there arises additional dependence of the conditions for the stability on the step size of discretization.

3.   A non-standard discretized model

Let us introduce the discrete version of system (2.2) with the use of the NSDM, following the approach presented in ^[3].

In our case applying the approach from ^[3] means replacing a bilinear term $x_n y_n$ occurring in system (2.4) with a term $x_{n+1} y_n$ or $x_n y_{n+1}$ in the first, second or both equations. In order to keep $x_n$ terms (meaning the scaled number of healthy people) positive, we introduce $x_{n+1} y_n$ in the first equation. The positivity of $y_n$ terms is easily preserved without any modifications of $x_n y_n$ in the second equation. Based on these assumptions, we modify system (2.4) and obtain the system having form

Transforming Eq (3.1a), we write system (3.1) as

Now we have the clear dependence of the variables' values for the $(n+1)$-th step on the variables' values for the previous one.

3.1.   Basic properties

Let us first present some basic properties of system (3.2).

Basing again on the approach from ^[20], we compute $\mathcal{R}_0$ for system (3.2) and obtained (2.3).

System (3.2) is a discrete dynamical system described by the function

that is the orbit of any arbitrary point $(x_0, y_0)$, indicated by $\mathcal{O}$, has the form $\mathcal{O} = ((x_0, y_0), H (x_0, y_0),$ $H^2 (x_0, y_0), \ldots).$ Because of the meaning of the system variables, we should have $F^n(x, y) \geq 0$ and $G^n(x, y) \geq 0$ for nonnegative $x$ and $y$. Note that $F(x, y) > 0$ for inequality (Eq 2.6) held and $G(x, y) > 0$ for

Since $k > \mu$, inequality (Eq 3.4) is stronger than inequality (Eq 2.6). From system (3.3) we conclude the following

Corollary 1. Let inequality (Eq 3.4) hold and $x_0 \geq 0$. Then

1. if $y_0 > 0$, then $x_n, y_n > 0$ for $n \in \mathbb{N}_+$,

2. if $y_0 = 0$, then $x_n > 0$ and $y_n = 0$ for $n \in \mathbb{N}_+$.

One should note that thanks to inequality (Eq 3.4) and the non-negative initial sizes of both healthy and infected subpopulations, our model is biologically reasonable. The number of the healthy people is always positive and the number of the infected ones is always non-negative – it is equal to zero, however, in case of the lack of the infection in the population.

In the following we assume that inequality (Eq 3.4) holds and consider system (3.2) for $x, y \geq 0$.

Now we estimate values of the $x$ variable. This variable increases if $\tfrac{F(x, y)}{x} > 1$ for $x \ne 0$, that is equivalet to

giving

In the considered domain the function $\tilde{x}(y)$ has its maximal value for $y = 0$ equal to

Hence, in the following we assume $x \leq \frac{C}{\mu}$. Substituting this inequality in $F(x, y)$ yields

Since $\frac{C}{\mu} \gg 1$, $F(x, y)$ has the maximal value for $y = 0$ equal to

and we conclude:

Corollary 2. Let inequality (Eq 3.4) hold. If $x_q \leq \frac{C}{\mu}$ for a given $q \in \mathbb{N}$, then $x_{q+s} \leq \frac{C}{\mu}$ for any $s \in \mathbb{N}$.

In other words, if the number of healthy people in any iteration of system (3.2) is smaller than $\frac{C}{\mu}$, then this value will be an upper bound for further iterations. Having given $x_0 \leq \frac{C}{\mu}$, we rephrase Corollary 2.

Corollary 3. If inequality (Eq 3.4) holds and $x_0 \leq \tfrac{C}{\mu}$, then $x_n \leq \tfrac{C}{\mu}$ for any $n \in \mathbb{N}_+$.

In the following we consider system (3.2) in the invariant subset

Now we focus on the $y$ variable. For $y \ne 0$, this variable decreases if

giving

Note that if $\frac{C}{\mu} < k$, then inequality (Eq 3.7) is satisfied in the invariant subset $\Omega$.

Combining this analysis with the notion of $\mathcal{R}_0$, we get

Corollary 4. If $x_0 \leq \frac{C}{\mu}$ and $\mathcal{R}_0 < 1$, then for system (3.2) we have $y_{n+1} < y_n$ for $n \in \mathbb{N}$.

Now let us check if the size of the whole population is bounded from above. We rewrite system (3.2) to get difference quotients

Adding by sides the equations from system (3.8), we get

giving

We consider two cases:

$\bullet$ Let $x_{n+1} \geq x_{n}$. From Eq (3.10) we have

that is equivalent to

what appears in system (2.5). Analogously like in system (2.4) solving inequality (Eq 3.12) with inequality (Eq 2.6) gives inequality (Eq 2.7), what implies inequality (Eq 2.8). Naturally, if we additionally assume $w_0 \leq \frac{C}{\mu}$, we get $w_n \leq \frac{C}{\mu}$ for every $n \in \mathbb{N}_+$.

$\bullet$ Now let $x_{n+1} < x_{n}$ be true.

First we assume that $\mathcal{R}_0 < 1$ and $x_0 \leq \tfrac{C}{\mu}$. From Corollary 4 we have $y_{n+1} < y_n$, and with $x_{n+1} < x_{n}$ it gives $w_{n+1} < w_n$. Hence the population size is bounded from above.

If $\mathcal{R}_0 > 1$ and $x_n < k$, then $y_{n+1} < y_n$, giving again the boundedness of the population size from above.

The case $\mathcal{R}_0 > 1$ and $x_n > k$ is complicated and we will omit its analysis in this paper.

Graphs in Figure 1 present the boundedness of the whole population which dynamics is described by system (3.2). The initial condition $(x_0, y_0) = (0.7; 0.2)$ and the values of parameters $\mu = 0.1$, $\alpha = 0.4$ and $h = 0.1$ were chosen. For Figure 1a we assumed that $C = 0.08$ and for Figure 1b we chose $C = 0.25$ so that $\mathcal{R}_0 < 1$ and $\mathcal{R}_0 > 1$ are fulfilled, accordingly. For each inequality $300$ and $600$ first iterations of system (3.2) were conducted, respectively.

3.2.   Stability of stationary states

Let us focus on stability analysis. System (3.2) has the same stationary states, and consequently the same conditions for their existence, as in systems (2.2) and (2.4). Jacobian matrix of system (3.2) has the following form

Let us start from the stability of $E_d$.

Theorem 1. If inequality (Eq 3.4) holds, then the disease-free state $E_d$ of system (3.2) is a sink for $\mathcal{R}_0 < 1$, non-hyperbolic for $\mathcal{R}_0 = 1$ and a saddle for $\mathcal{R}_0 > 1$. Moreover, if $x_0 < \frac{C}{\mu}$ and $\mathcal{R}_0 < 1$, then $E_d$ is globally stable.

Proof. Inequality (Eq 3.4) guarantees the positivity of $x$ and $y$. The Jacobian matrix for the $E_d$ state has the form

The conditions for $E_d$ stability are $|\lambda_i| < 1$, $i = 1$, $2$. Note that $\lambda_1 = 1-h\mu$ satisfies it under the assumption inequality (Eq 3.4). On the other hand, $\lambda_2 = 1-hk+h \tfrac{C}{\mu}$ is positive for inequality (Eq 3.4) and $\lambda_2 < 1$ holds if $\tfrac{C}{\mu} < k$, that is $\mathcal{R}_0 < 1$.

Now we show the global stability of $E_d$ for $x_0 < \frac{C}{\mu}$ and $\mathcal{R}_0 < 1$. From Corollary 3 we have $x_n \leq \tfrac{C}{\mu}$. Moreover, Corollary 3 and Eq (3.2b) give

From $\mathcal{R}_0 < 1$ and inequality (Eq 3.4) we get $0 < 1-h k \left(1 - \mathcal{R}_0 \right) < 1$ and

Now let us add by sides the equations from system (3.1). We get

This yields

Let us take $\varepsilon > 0$. Since Eq (3.14) holds, we state that

Combining inequalities (Eq 3.16) and (Eq 3.17), we get for $n > N_1$

giving

Becuase

we can choose sufficiently large $n$ ($> N_1$) and we get

Using Eq (3.15) again we obtain

From Eq (3.14) for sufficiently large $n$ we have

and

yielding

We state that for sufficiently large $n$

holds. Combining inequalities (Eq 3.18) and (Eq 3.19), for sufficiently large $n$ we get

implying

Clearly, we have

Both Eqs (3.14) and (3.20) yield

proving global stability of $E_d$.

Note that conditions of stability of $E_d$ do not depend on the step size of discretization $h$, while this is not the case of the endemic state $E_e$ we discuss below.

Figure 2 includes a phase portrait illustrating the global stability of $E_d$ for system (3.2). The values of parameters $C = 0.05$, $\mu = 0.1$, $\alpha = 0.4$ and $h = 0.1$ were chosen so that the condition $\mathcal{R}_0 < 1$ is fulfilled. Four initial conditions $(x_0, y_0) = (0.8;0.2)$, $(x_0, y_0) = (0.15;0.05)$, $(x_0, y_0) = (0.2;0.25)$ and $(x_0, y_0) = (0.3;0.5)$ were chosen. For each initial condition $50000$ iterations were conducted.

Now we investigate the stability of the endemic stationary state $E_e$. Hence, in the following analysis we assume $\mathcal{R}_0 > 1$. Let us define parameters

and

Furthermore, we will refer to the $\delta$ value presented in Eq (2.12). See that because $C \gg \mu$, we have $h_1 > 0$ for $\mathcal{R}_0 > 1$. Positivity of $h_2$ is guaranteed for $\delta > 0$. In order to present the results transparently, we first consider the case $\delta > 0$. The condition $\delta < 0$ will be discussed later. The case $\delta = 0$ is not generic and will be not analyzed.

Let us formulate the theorem.

Theorem 2. Suppose that the state $E_e$ of system (3.2) exists. Moreover, assume that $\delta > 0$. The state $E_e$ is then a sink for $h < h_2$, a spiral sink for $h \in \left(h_2, h_1 \right)$ and $h_2 < h_1$, a spiral source for $h > h_1 > h_2$, and non-hyperbolic for $h = h_1 > h_2$, where $h_1$ and $h_2$ are defined in Eqs (3.22) and (3.23), respectively.

Proof. The Jacobian matrix of system (3.2) for $E_e$ has the form

To simplify the notation, we define

From inequality (Eq 2.6) and $\mathcal{R}_0 > 1$ we have $\upsilon > 0$ and $\varpi > 0$. Using Eqs (3.24) and (3.25), we rewrite $M(E_e)$ as

The characteristic polynomial of $M(E_e)$ reads

while eigenvalues are equal to

where

Let us check for which $h$ a condition $\widehat{\delta} > 0$ holds. We get the sequence of equivalent inequalities

leading to

From $\delta > 0$ we have the positivity of the right-hand side of inequality (Eq 3.28). Combining inequality (Eq 3.28) and Eq (3.23), we get $h < h_2$.

Analogously we have $\widehat{\delta} < 0$ for

The further analysis considers different signs of $\widehat{\delta}$. We omit the case $\widehat{\delta} = 0$.

$\bullet$ Let us consider $\widehat{\delta} > 0$, for which we get the real eigenvalues such that $\lambda_1 > \lambda_2$. Note that Polynomial Eq (3.27) has two real zeros and the coefficient from its linear term is negative. Hence, $\lambda_1, \lambda_2 > 0$ and for the local stability of $E_e$ it is enough to fulfill $\lambda_1 < 1$. Looking for this inequality we obtain

giving

The right-hand side of inequality (Eq 3.30) can be written as

The above expression is always positive, so we can square both sides of inequality (Eq 3.30) obtaining

which is obvious from the definition of $\widehat{\delta}$.

From $\lambda_1 < 1$ we do not get additional conditions for the local stability. Hence, we state that existing $E_e$ is locally stable if $\delta > 0$ and $h < h_2$ hold.

$\bullet$ Now we consider the case $\widehat{\delta} < 0$, making the eigenvalues complex. We have

The condition $|\lambda| < 1$ can be written as

giving

We denote the left-hand side of inequality (Eq 3.32) as $\bar{P}(h)$, which is a quadratic trinomial with a negative constant term. From the existence of $E_e$ we have $y_e (C - \mu k) > 0$. Hence, $\bar{P} (h)$ has one positive zero equal to

From Eq (3.22) we state that this zero is equal to $h_1$. We conclude that $h < h_1$ should hold in order to fulfil $\bar{P} (h) < 0$. Combining $h < h_1$ and inequality (Eq 3.29), we get $h \in \left(h_2, h_1 \right)$, being true for $h_2 < h_1$.

Let us focus again on Theorem 2. Note that after exceeding $h_2$ the state $E_e$ changes from a sink to a spiral sink, while the local stability is preserved. Moreover, $h_1$ is a critical value in the context of stability – after exceeding $h_1$ the state $E_e$ stops being stable.

Now let us assume $\delta < 0$. Then we have $h_2 < 0$, corresponding to $\widehat{\delta} < 0$. Using that and Theorem 2, we conclude:

Corollary 5. Suppose that the stationary state $E_e$ of system (3.2) exists. Let us also assume that $\delta < 0$. The state $E_e$ is then a spiral sink for $h < h_1$, a spiral source for $h > h_1$ and non-hyperbolic for $h = h_1$.

3.3.   Version of Theorem 2 and Corollary 5 taking into account positivity of solutions.

Now we assume that inequality (Eq 3.4) is true, guaranteeing the positivity of solutions of system (3.2). We check conditions for which the threshold values of $h$ defined in the previous subsection satisfy this inequality. Let us first check the condition giving

After some calculations this inequality is equivalent to

where

Note that $a_0 > 0$. Let $\delta_C$ be a discriminant of $W(C)$. After some tedious calculations we obtain

Because $\delta_C > 0$, we have two real zeros of $W(C)$ equal to

The inequalities $a_2 > 0$ and $a_1 < 0$ can be written as

and

Let us take $m \in \left(0, \frac{1}{3} \right)$. Assuming inequality (Eq 3.39), from inequality (Eq 3.40) we have

We have $1 - 2 m > 0$, so the last inequality from above is always true.

Hence, $a_2 > 0$ implies $a_1 < 0$. We also see that $a_2 > 0$ and $a_1 > 0$ cannot hold simultaneously.

The fulfillment of inequality (Eq 3.35) depends on the sign of $a_2$ solely. Because inequalities (Eq 3.34) and (Eq 3.35) are equivalent, we state that:

1) for $a_2 > 0$ inequality (Eq 3.34) holds for $C > 0$ and $C \in (C_a, C_b)$,

2) for $a_2 < 0$ inequality (Eq 3.34) holds for $C > C_b > 0$.

Now let us focus on the inequality

Using the $y_e$ notation, we rewrite $h_1$ as

Substituting Eq (3.42) into inequality (Eq 3.41), we get

that is

Inequality (Eq 3.43) is justified for

If the condition contrary to Eq (3.44) holds, then we have

Let us assume that Eq (3.44) is true. We can then square both sides of inequality (Eq 3.43) obtaining

what gives

Note that from $C \gg \mu$ we have $k (1 + \mu) (C-\mu) \gg 0$. Hence we state that inequality (Eq 3.46) never holds. It means that inequality (Eq 3.41) is contradictory and we shall assume inequality (Eq 3.45).

Now let us consider Theorem 2 and Corollary 5 in the context of positivity of solutions of system (3.2). Let us emphasize again that the condition $h > h_1$ does not hold. We consider the sign of $a_2$ defined in Eq (3.36). From Theorem 2 we obtain the following corollary:

Corollary 6. Assume that the stationary state $E_e$ of system (3.2) exists and define $C_a$ and $C_b$ as in Eq (3.38). If $\delta > 0$ and one set of the conditions:

or

holds, then $E_e$ is a sink for $h < h_2$ and a spiral sink for $h \in \left(h_2, \tfrac{1}{k} \right)$. Moreover, the solutions of system (3.2) are positive for $x_0 \geq 0$ and $y_0 > 0$.

Let us look at Corollary 5. Remind that $\delta < 0$ is equivalent to $h_2 < 0$. We see that then inequality (Eq 3.34) is always true and we should only consider inequality (Eq 3.45). We reject the case when $h > h_1$. Hence we have:

Corollary 7. Let the state $E_e$ of system (3.2) exist and inequality (Eq 3.4) be true. Then $E_e$ is a spiral sink and the solutions of system (3.2) are positive for $x_0 \geq 0$ and $y_0 > 0$.

3.4.   Possibility of Neimark-Sacker bifurcation appearance

Let us consider the possibility of the Neimark-Sacker bifurcation appearance for the $E_e$ state in system (3.2). The necessary condition giving this possibility is $|\lambda_1| = |\lambda_2| = 1$. Using the calculations presented above, we determine the necessary conditions for NSB:

Corollary 8. Let us assume that either $\delta > 0$ and $\widehat{\delta} < 0$ or $\delta < 0$. Neimark-Sacker bifurcation for the state $E_e$ in system (3.2) can happen if $h = h_1$.

Determining the conditions for NSB in systems based on the NSDM is in general difficult. Because of that, we only numerically illustrate the possibility of the NSB appearance for system (3.2).

We conduct numerical simulations which results are presented in Figures 3 and 4. In each figure first $50,000$ iterations of system (3.2) were presented. For simulations the values of parameters: $C = 0.19$, $\mu = 0.1$ and $\alpha = 0.3$ were taken. The values were chosen so that the condition $\delta < 0$ is fulfilled. For Figure 3 $(x_0, y_0) = (1.52; 0.2)$ and $h = 4.5$ were taken. In case of Figure 4 the step size $h = 4.65$ was assumed and for each Figures 4a and 4b the initial condition $(x_0, y_0) = (1.52; 0.2)$ and $(x_0, y_0) = (1.49; 0.3)$ was respectively chosen. Note that for $h = 4.5$ the state $E_e$ is a sink and for $h = 4.65$ we get the invariant curve attracting the solutions from both inside and outside the curve, what is shown in Figures 4a, b. Analysing the behavior illustrated in Figures 3 and 4, we guess that in system (3.2) there is a possibility of NSB for the endemic stationary state for $h = h_1$.

4.   Conclusions

In this paper we analyzed the non-standard discretized system describing epidemic spread in the population consisting of two groups: susceptible (healthy) and infected people. We presented basic properties of the assumed system of difference equations, focusing on positivity and boundedness of the variables. In both discrete model (2.4) with the EEM and system (3.2) with the NSDM we obtained boundedness of the variables with a proper initial condition (besides the unchecked case $\mathcal{R}_0 > 1$ and $x_n > k$ for the NSDM). We also proved the global stability of the disease-free stationary state $E_d$ and the local stability of the endemic state $E_e$.

The results for $E_d$ for sufficiently small step of discretization are analogous to those for this state in case of the continuous model. Note that these results are nearly the same as for system (2.4).

The state $E_e$ for the continuous case is globally stable if it exists (see ^[17]). In case of the NSDM we proved the local stability of $E_e$ for $h < h_1$. These results are in accordance with the fact that the behavior of a discrete system for a sufficiently small step of discretization is analogous to the behavior of its continuous counterpart. Even the analysis of local stability of the endemic state occurred complex for system (3.2), so we expect that global stability analysis would be much more complex, which is in general the case for systems based on the NSDM. This is mainly related to the dependence of stability of $E_e$ on the step of discretization and the possibility of bifurcations with respect to this parameter. Note that that behavior of systems (3.2) and (2.4) is analogous near accordingly $h_1$ and $h_c$. We conclude that in both discrete models for the endemic state $E_e$ there is a critical value of the step size of discretization and after exceeding this value the state $E_e$ loses the local stability. After necessary computations, we got that $h_1 > h_c$. It means that the loss of $E_e$ stability happens in system (3.2) for a bigger value of $h$ than in system (2.4). Hence, we state that using the NSDM in a discrete system gives better approximation of the analogical continuous counterpart comparing to using the EEM.

Comparing the results for the continuous model to those for both discrete ones, we conclude that the conditions for the stability of the stationary states requires the bounded value of the step size of discretization from above for the discrete models. In case of the NSDM the conditions are weaker than in case of the EEM.

We checked the possibility of the Neimark-Sacker bifurcation appearance in the analyzed system. We presented one specific example for which an isolated invariant closed curve seems to appear. Looking at Figures 3 and 4, we state that the critical value for the NSB belongs to the interval $(4.5; 4.65)$. Note that after exceeding the critical value, the stationary state loses its stability and there arises a stable closed invariant curve surrounding this state. We conclude that in system (3.2) there can be a supercritcal NSB. We conducted numerical simulations for possibility of NSB appearance in the discrete model with the EEM. The simulations suggest the appearance of a subcritical NSB for the same stationary state $E_e$, differently to the model with the non-standard discretization. In order to prove that for two types of discretization we have different types of NSB, mathematical analysis should be conducted.

Acknowledgments

We would like to thank the constructive feedback provided by the reviewers.

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this article.