DYNAMICS OF AN ULTRA-DISCRETE SIR EPIDEMIC MODEL WITH TIME DELAY

. We propose an ultra-discretization for an SIR epidemic model with time delay. It is proven that the ultra-discrete model has a threshold property concerning global attractivity of equilibria as shown in diﬀerential and diﬀerence equation models. We also study an interesting convergence pattern of the solution, which is illustrated in a two-dimensional lattice.

1. Introduction. Differential equations are used for modeling phenomena in many disciplines e.g. economy, physics, engineering, biology and chemistry. Since exact (analytical) solutions are not available in many cases if the mathematical model is given as a nonlinear system, numerical solutions have enhanced our understanding of the mathematical model [7]. Traditional numerical schemes such as Euler's method and Runge-Kutta method, however, may induce numerical instability, as the discretization would change properties of solutions, such as stability and positivity, of the original model, thus numerical scheme should be carefully chosen to preserve nature of the original system, see e.g. [4,8].
The authors in [10] propose a system of difference equations as a discrete counterpart of a continuous epidemic model, describing disease transmission dynamics in continuous time. Qualitative properties of the model such as positivity, boundedness and convergence of the solutions are investigated. As a continuation work of [10] in [11] a general system of difference equations is analyzed. In both papers, proving convergence of the solution, when a parameter called the basic reproduction number is greater than one, seems to be a challenging problem, while it is known that the corresponding continuous model has an equilibrium that is globally stable. A "good" discretization is found in [23,24,5] for a class of epidemic models formulated by delay differential equations in [2,3,29]. The authors in [24] prove uniform persistence of the solution, corresponding to a result in [29] for a continuous SIRS model. The authors in [5] prove global stability of the endemic equilibrium by a Lyapunov function, corresponding to the result established in [14] (see also [6]). The discretization used in [24,5] is a variation of backward Euler's method and is indeed known as Micken's nonstandard finite difference method [15]. See also [12,16] for the application of Micken's nonstandard finite difference method to ordinary differential equation models. It is also known that some discrete-time epidemic models exhibit periodic and chaotic behavior [1].
Ultra-discretization is proposed as a procedure to obtain a discrete system, where unknown variables also take discrete values, thus a cellular automaton is defined, see [21,19]. In [27,22,20] the authors study discrete and ultra-discrete models for epidemic models. In those papers it is shown that simple ultra-discrete models can capture the disease transmission dynamics, which is seen in the original differential equation models. In [27] the authors find conserved quantities for some special cases. Cellular automata have been used to model disease transmission dynamics, see e.g. [25,26] and references therein. Since cellular automata are computational models, in general, it is not straightforward to perform a mathematical analysis, in order to provide theoretical insights into simulation studies. Our analytical approach in this paper for the ultra-discrete model may be used to complement numerical simulation studies for some cellular automaton models.
In this paper we start with a special case of the model considered in [5,24] for a discrete analogue of an epidemic model formulated by a system of delay differential equations in [2,3,29,14]. As delay differential equation form an infinite dimensional dynamical system [9], after discretizing the system, we obtain a system of difference equations of higher order, which is slightly complicated compared to the model considered in [27]. From such a model we derive an ultra-discrete model, which is formulated as a couple of piecewise linear difference equations (see e.g. [13]). Due to time delay, integrability can not be usually expected and the application of the ultra-discretization to a non-integrable system is still limited, but see also [17,18]. We here prove that the ultra-discrete model exactly has the same threshold property regarding global attractivity as in [5,14]. In Section 4 we visualize the convergence of the solution in a two-dimensional lattice and observe some interesting convergence patterns. For a special initial condition, we derive a simple recurrence relation for the solution in Lemma 3. The relation can explain the illustrated convergence pattern.
The paper is organized as follows. In Section 2 we introduce a system of difference equations, which is a special case of the model studied in [5,24] for discrete analogue of continuous models. Applying a variable transformation together with taking a limit, we derive an ultra-discrete model. In Section 3 global behavior of the solutions is discussed. We prove that the model exhibits the threshold behavior, similar to the difference equation studied in [5]. We here find that a subsequence of the solution has a monotone property and this monotonicity is used for the proof. In Section 4 we illustrate the solution behavior in a two-dimensional lattice. To explain the convergence pattern, we consider a special initial condition and derive a simple recurrence relation. We summarize our results in Section 5.
2. Ultra discretization of an epidemic model. Our starting point is the following system of difference equations with a positive initial condition. System (1) is a special case of the model considered in [24,5]. See also [6] for a model with nonlinear incidence. Here M, B and Γ are positive parameters and ω is a positive integer. Let us define One can prove that system (1) always has the disease free equilibrium given by (1, 0, 0) and that there exists an endemic equilibrium given by The authors in [5] prove the following threshold type behavior.
(i) If R 0 ≤ 1 holds, then the disease free equilibrium of (1) is globally asymptotically stable. (ii) If R 0 > 1 holds, then the disease free equilibrium of (1) is unstable and the endemic equilibrium is globally asymptotically stable.
System (1) is proposed as a discrete analogue of the following disease transmission dynamics model in continuous time: where S(t), I(t) and R(t) respectively denote fraction of susceptible, infective and recovered population at time t. The constant M > 0 denotes the death rate. The constant B > 0 is transmission coefficient and the constant Γ > 0 is the recovery rate. The non-negative constant τ ≥ 0 can be interpreted as incubation period of infection. This model was developed in [2] to describe transmission dynamics of a vector-borne disease mosquito, see also [3]. In [14] it is shown that the continuous model (2) exhibits the same threshold behavior as in Theorem 2.1 for (1). Let us now derive an ultra-discrete model from system (1), following the same procedure as in [27,22,20]. Since the third equation of (1) does not appear in the first and second equations for S and I of (1), we focus on the first and second equations of (1a) and (1b). For > 0 we introduce variables x and y via S n = e xn/ and I n = e yn/ , and parameters µ, β and γ through M = e µ/ , B = e β/ and Γ = e γ/ .
Notice that (1a) and (1b) are equivalently written as Applying the variable transformation to (3) with letting → +0, we get the following ultra-discrete model The key relation used here is the following limit The initial condition of the system (4) is given as Moreover, parameters µ, β and γ are given as real numbers. For simplicity, we assume that 0 < µ ≤ γ. From (4) we then get 3. Global attractivity of equilibria. In this section, we study the global asymptotic behavior of the solutions of (6).
Lemma 3.1. For any solutions, there existsn ∈ N + such that x n ≤ 0 and y n ≤ µ−γ for n ≥n.
Proof. Let us assume that x n ≥ µ for some n ≥ 0. Then from (6a) one has that thus x n is decreasing with respect to n as long as x n ≥ µ. Therefore, there exists k such that x k−1 ≥ µ and x k < µ. Then from (6a), it follows Thus we get that x m ≤ 0 for all m ≥ k + 1. So hereafter, without loss of generality we can set x n ≤ 0 for all n ≥ 0. Since we have follows. Therefore, from (6b), we have If y n ≤ µ then it is easy to see y m ≤ µ − γ for all m ≥ n + 1. Let us assume that y n ≥ µ for some n ≥ 0. Then we have that y n+1 ≤ y n − γ, thus y n is decreasing with respect to n as long as y n ≥ µ. Therefore, there exists such that y −1 ≥ µ and y < µ. Then one obtains y +1 ≤ µ − γ. According to the previous discussion, there existsn ≥ + 1 such that y m ≤ µ − γ for all m ≥n.
From Lemma 3.1, without loss of generality, we can set the initial condition as Note that Lemma 3.1 implies is an invariant set. To discuss global attractivity of equilibria, we now consider the following system in (8). Proof. Since for any n ≥ 0 one has that β + y n−ω ≤ β + µ − γ < µ from Lemma 3.1, it follows max (µ, β + y n−ω ) = µ.
We now show that every solution converges to the non-trivial equilibrium. Proof. Let y := y (ω+1) , y (ω+1)−1 , · · · , y (ω+1)−ω for ∈ N + . From Proposition 2 one can see that i.e. each component of y converges to the equilibrium as → ∞. Then, there exists a sufficiently large integer m such that y n = · · · = y n−ω = µ − γ for n ≥ m. For n ≥ m we obtain Corresponding to the first and second parts of Theorem 2.1, we show the threshold behavior in Theorems 3.2 and 3.3 for the ultra-discrete epidemic model (9).
In Figures 1 and 2 we plot x n and y n with respect to n. In Figure 1, we set ω = 0 so that the ultra-discrete model (9) has no time delay. The initial condition is given as x 0 = −3, y 0 = −13. We set the parameters as µ = 1, γ = 6 and β = 3 in Figure 1(a) while µ = 1, γ = 6 and β = 9 in Figure 1(b). As in Theorems 3.2 and 3.3, one can see that y tends to −∞ as n → ∞ in Figure 1(a) and that y tends to µ − γ as n → ∞ in Figure 1     In Figure 2, we set ω = 10. The initial condition is chosen as x 0 = −3 and y −j = −16 + j for j ∈ {0, 1, · · · , ω}. We set the same parameters as Figures 1(a)  and 1(b) in Figures 2(a) and 2(b), respectively. As in Theorems 3.2 and 3.3, again one can see that y tends to −∞ as n → ∞ in Figure 2(a) and that y tends to µ − γ in Figure 2(b). Comparing Figure 1 with Figure 2, it can be seen that the solution is monotone for ω = 0. In the ultra discrete model (9) time delay does not change qualitative dynamics, but changes the solution behavior. 4. Monotone convergence in a two-dimensional lattice. We here visualize the convergence of a solution using a two-dimensional lattice. We consider the case that β − γ > 0 holds, so the solution converges to the non-trivial equilibrium (−β + γ, µ − γ), according to Theorem 3.3. The two-dimensional lattice is constructed using two variables: time step and time delay. More specifically, let u −j m := y m(ω+1)−j for m ∈ N + and j ∈ {0, 1, · · · , ω} to set a sequence with two variables (j, m). From Proposition 2, we have for j ∈ {0, 1, · · · , ω − 1} and In Figure 3, specifying parameters and the initial condition, we write a numerical value of u −j m in the corresponding lattice. The black lattice represents that u −j m = µ − γ i.e., the solution reaches the equilibrium. Of course, from Theorem 3.3, the solution converges to the homogeneous equilibrium (with respect to j) as m → ∞. Moreover, one may observe interesting stepwise shape representing the convergence pattern of the solution.
To explain the pattern we consider a special solution and derive an explicit recurrence relation for the solution in Proposition 3. For k ∈ {1, 2, · · · ω − 1} set the initial condition as with p < q < µ − β.
When m = 0 the k-th component is the biggest component, i.e., u −k 0 = max 0≤j≤ω u −j 0 . In Lemma 4.1 we show that this relation is preserved when m = 1. It is also shown that , 2, · · · , k} . In this way u −j 1 is ordered for j ∈ {0, 1, · · · , k}. Furthermore, from Lemma 3, one can see that the monotone ordering is preserved for m ≥ 1. Consequently, the stepwise shape appears as the convergence pattern illustrated in Although here we consider a special solution such that the initial condition is given as in (14), we can show that linear combination of two solutions can produce a complicated pattern for the convergence of the solution. Consider two solutions u −j m and v −j m of (12) and (13). Then is shown to be a solution of (12) and (13) as follows. Assume that u −j m−1 and v −j m−1 are less than µ − β, then w −j m−1 < µ − β. Using (12) one can compute Assume that either u −j m−1 > µ − β or v −j m−1 > µ − β. Then w −j m−1 > µ − β. In this case one can prove w −j m = µ − γ. Let the solution depicted in Figure 3(a) be u and the solution depicted in Figure  3(b) be v. Convergence pattern of the solution w, which is given by (27), is illustrated in Figure 3(c). From our discussion, the convergence pattern in Figure 3(c) appears as the composition of the stepwise convergence pattern which is observed in Figures 3(a) and 3(b). 5. Conclusion. In this paper we consider an ultra-discrete model with time delay. The model is derived from a discrete epidemic model studied in [5,24]. In Theorems 3.2 and 3.3, we show that the ultra-discrete model also has the threshold property concerning global attractivity of equilibria, similar to the discrete epidemic model [5] and the continuous epidemic model [14]. For the proof of global attractivity of the non-trivial equilibrium in Theorem 3.3, we reduce the system (9) to the scalar difference equation in Proposition 2 and then use a certain monotone property of the solution, which is our important finding.
In Section 4 we further derive a simple recurrence relation for the solution, assuming a special condition for the initial condition. The relation derived in Proposition 3 clearly shows that each component monotonically increases towards the equilibrium. Two-dimensional lattice seems to be an informative tool to illustrate such a convergence patten.
For the ultra-discrete model (9), a linear combination of two solutions is shown to be a solution. It would be interesting to investigate the structure of the solutions of ultra-discrete models. Figure 3 may remind of an elementary cellular automaton (see e.g. rule 252 in [28]). Exploring a possible connection to elementary cellular automaton is our future work.