Fractional-Order Discrete-Time SIR Epidemic Model with Vaccination: Chaos and Complexity

: This research presents a new fractional-order discrete-time susceptible-infected-recovered (SIR) epidemic model with vaccination. The dynamical behavior of the suggested model is examined analytically and numerically. Through using phase attractors, bifurcation diagrams, maximum Lyapunov exponent and the 0 − 1 test, it is veriﬁed that the newly introduced fractional discrete SIR epidemic model vaccination with both commensurate and incommensurate fractional orders has chaotic behavior. The discrete fractional model gives more complex dynamics for incommensurate fractional orders compared to commensurate fractional orders. The reasonable range of commensurate fractional orders is between γ = 0.8712 and γ = 1, while the reasonable range of incommensurate fractional orders is between γ 2 = 0.77 and γ 2 = 1. Furthermore, the complexity analysis is performed using approximate entropy ( ApEn ) and C 0 complexity to conﬁrm the existence of chaos. Finally, simulations were carried out on MATLAB to verify the efﬁcacy of the given ﬁndings.


Introduction
Epidemiology is a topic of research in the biological sciences, which explores all the elements that influence whether there are diseases and disorders. Over the last few years, the outbreaks of several diseases such as SARS, Ebola, and COVID-19 have increased the attention of many researchers to the study of epidemiology. In this regard, the construction of infectious diseases mathematical models and discussion of their dynamics are very significant, since these models are used to aid prevent and control diseases. Numerous studies on mathematical epidemic models have been examined for a variety of diseases. These models may be classified as either continuous-time [1][2][3][4][5] or discrete-time models. Differential equations are used to create continuous-time models, while difference equations are used to formulate discrete-time models. Recently, discrete models have been extensively employed to evaluate and study infectious diseases as opposed to continuous models since epidemic data are available during discrete time periods. Furthermore, the discrete-time epidemic models display more complicated dynamics than the continuous-time epidemic models. Numerous significant types of discrete models may be found in [6][7][8][9].
Vaccination is a significant component of disease prevention strategies across the globe. A vaccine's effectiveness is based on its ability to reduce the number of susceptible individuals in order to prevent infectious diseases from spreading through the population [10]. The effect of vaccination on the dynamics of an epidemic model has been studied by adding a component for the vaccinated individuals to the epidemic model and obtaining the epidemic model with vaccination. For example, in [11] the authors studied the dynamical behavior of a discrete SIR epidemic model with a constant vaccination strategy. The stability and bifurcations of a discrete susceptible-infected-susceptible (SIS) epidemic model with vaccination were investigated in [12], while Xiang et al. in [13] developed a discrete SIRS model that included vaccination and examined its dynamical behavior.
Fractional calculus is a subject in mathematical analysis, where it can be considered as a generalization of integer calculus [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. Despite this, only in the past decades has it been extensively examined, owing to its broad range of use in many areas. Fractional-order derivatives are relatively more accurate than integer-order derivatives because they serve as an effective tool for describing the memory effect throughout all types of processing and materials [29]. Memory effects refer to the fact that the states of systems with a fractional order are determined by all previous states. Over the last several years, numerous academics have concentrated their efforts on discrete fractional calculus [30][31][32]. Following the proposal of that special field, researchers have become more interested in its applications in neural networks, physics, biology, etc. Meanwhile, since Wu et al. [33] suggested the first chaotic fractional logistic map, the dynamics and control of various fractional-order chaotic discrete-time systems have been intensively studied [34][35][36][37][38].
Fractional-order derivatives have been extensively used in epidemiology, with the majority of these studies focusing on SIR-type models [39][40][41][42][43][44]. In [45], Selvam et al. investigated the stability of the fractional-order SIR epidemic model using a discretization process, whereas Naik in [46] studied the global dynamics of a fractional-order SIR epidemic model with memory. The analysis and numerical solution of a novel fractional-order SIR dengue model was investigated in [47]. In [48], the stability analysis and bifurcation control for a delayed fractional-order SIR epidemic model with incommensurate order was considered. It is worth noting that all of these studies are based on continuous-time systems. The motivation behind this work, particularly with the current spread of the COVID-19 pandemic, is to assess the benefits of the fractional discrete epidemic model in the field of epidemiology. So far, to the best of our knowledge, the dynamic analysis of a fractional-order discrete-time epidemic model with vaccination based on a Caputo-like difference operator has not yet been investigated. This piqued our interest and inspired us to study the phenomenon and investigate the behavior of a fractional SIR epidemic model with vaccination when the fractional orders are commensurate and incommensurate.
The goal of this article is to contribute to the field of epidemiology by introducing a novel discrete-time SIR epidemic model with vaccination with both commensurate and incommensurate fractional orders. The work is organized as follows. In Section 2, we present the integer-order form of the discrete SIR model and we recast the mathematical model to take the form of the fractional-order discrete SIR model. In Section 3, we look at how the range and type of chaotic behaviors are affected by commensurate and incommensurate fractional orders. Furthermore, we use the 0−1 test method to distinguish between regular and chaotic behavior. In Section 4, we use the C 0 complexity and approximate entropy (ApEn) to analyze the complexity of the proposed discrete SIR epidemic model. Throughout the paper, numerical simulations are used to demonstrate and verify the results.

Mathematical Model
In this section, we introduce the discrete SIR model as an integer-order discrete system, and then we reformulate it as a fractional-order discrete system by employing the Caputoleft difference operator.

Integer-Order Discrete Model
Consider the following description of the SIR epidemic model described in [49], in which the population densities S m , I m , and R m represent the number of susceptible, infected, and recovered individuals at time m, respectively: where α > 0 is the contact number, 0 < β < 1 is is the probability of birth, and 0 < σ < 1 is the probability of recovery. The flow diagram for model (3) is seen in Figure 1. In [50], the authors focused on the dynamics of an SIR epidemic model by incorporating vaccination in the model (1). The proportion of persons who had been vaccinated was represented by the parameter p, where 0 < p < 1; and the remainder had a risk of susceptibility to infection in the proportion 1 − p. The following is the general SIR epidemic model with vaccination, which takes the following form: (S 0 , I 0 , R 0 ) are initial conditions that are positive real values, with S 0 + I 0 + R 0 = N. N is the total population size. Furthermore, by employing the relation S m + I m + R m = N, we replace R m by N − S m − I m and we get the following system: Observing that the two equations concerning (S, I) in the SIR epidemic model (3) do not contain R, they are thus independent of the third one. As a result, the three-dimensional SIR epidemic model (3) can be reduced to a two-dimensional model, which results in the following equivalent system: when α = 4, β = 0.8, σ = 0.1, p = 0.005, and N = 100, the discrete SIR epidemic model with vaccination (4) exhibits a chaotic attractor as shown in Figure 2.

Fractional-Order Discrete Model
In this study, a discrete-time SIR epidemic model with vaccination with both commensurate and incommensurate orders is considered. Herein, the first-order difference of the SIR model (4) is formulated as: As mentioned above, we use the Caputo-left difference operator to create the fractional version of the discrete-time SIR epidemic model with vaccination. The Caputo-left difference operator is defined as follows: Definition 1. [51] The γ-Caputo fractional difference operator for a function h(s), is defined by and ∆ m h(τ) are the falling factorial function and the m-th integer difference operator, respectively, which are defined as and Remark 1. If we consider that m = 1, we can define the γ-Caputo fractional difference operator by Definition 2. [30] The γ-th fractional sum for a function h is defined as Now, the discrete SIR epidemic model with vaccination (4) may be expressed in fractional order as follows: 2,3) are the fractional order such that 0 < γ i < 1. Note that if γ 1 = γ 2 , then system (11) is referred to as a commensurate fractional-order system, whereas it is referred to as an incommensurate fractional-order system if γ 1 = γ 2 .

Dynamical Analysis and Numerical Simulations
In this section, we investigate whether the previously suggested fractional-order discrete-time SIR epidemic model with vaccination (11) exhibits chaotic behavior in both the commensurate and incommensurate fractional orders. This study is conducted utilizing many numerical techniques, including bifurcation diagrams, Lyapunov exponent computations, plotting of phase portraits in S-I projection and applying the 0−1 method.
In order to show these, we first present a theorem that enables us to obtain the numerical formula for the discrete fractional model that is subsequently discussed. Theorem 1. [52] For the fractional difference equation the unique solution of this initial value problem (1) is given by According to Theorem 1, the numerical formula of the fractional discrete-time SIR epidemic model (11) is designed as: Take a = 0 and since (s − τ + 1) (γ−1) is equal to Γ(s − τ)/Γ(s − τ − γ + 1), it follows from Theorem 1 that the numerical formula for (11) is obtained as: where S(0) and I(0) are the initial conditions. This is a novel type of fractional-order discrete-time SIR epidemic model with vaccination that possesses "memory effects". As seen in Equation (16), the states S(n) and I(n) are dependent on all previous variables S(0), S(1), . . . , S(n − 1), and I(0), I(1), . . . , I(n − 1).

Bifurcation Diagram and Maximum LEs
To investigate the dynamic of the commensurate fractional-order discrete-time SIR epidemic model with vaccination given in (11) regarding the fractional order γ = γ 1 = γ 2 , we set the parameters N = 100, β = 0.8, σ = 0.1, p = 0.005, and initial conditions (S(0), I(0)) = (70, 30). Figure 3 displays the phase space of model (11) for γ = 1. Figure 4 shows the states of the system for 150 points. Take note that when γ = 1, the fractional discrete model (11) refers to the integer-order model. Now, using the same initial condition and the same system parameters, the fractional discrete model for different fractional orders γ is presented in Figure 5. We see that as γ decreases, the phase portrait changes its shape between chaotic and periodic trajectories, whereas if we continue to decrease γ, the states of the fractional discrete model diverge to infinity.
The bifurcation diagrams of the fractional discrete model (11) where α varies in the interval [3, 4.5] is presented in Figure 6. Obviously, decreasing the fractional order has an effect on the interval in which chaos occurs, and the bifurcation diagram progressively shifts to the right. For example, when γ = 0.98, chaos begins to appear at α = 3.915 and the maximum chaotic range is reached at α = 4.128, whereas when γ = 0.91, chaos begins to appear at α = 3.849 and the maximum chaotic range is reached at α = 4.047.
To further study the influence of the commensurate fractional order on the dynamics of the fractional discrete-time SIR epidemic model (11), we plotted the bifurcation diagram for γ as a critical parameter and the findings are presented in Figure 6a. We notice that the fractional discrete SIR epidemic model (11) contains chaos and that the fractional order γ has an influence on the system's dynamics. When γ < 0.8712, the system (11) exhibits a transient state, which means that until a minimum number of iterations is reached, the proposed system's states approach a limited attractor before diverging to infinity.    An important tool used in fractional discrete systems is the Lyapunov exponents, which are used in conjunction with bifurcation diagrams to show chaos. As shown in [52], the Jacobian matrix algorithm is used to compute or estimate the maximum Lyapunov exponent, which is calculated in a similar way to the states in the fractional-order discrete model (11). The tangent map J i is defined as follows: Then, the Lyapunov exponents can be given by : where λ k are the eigenvalues of the Jacobian matrix J i . The method described above can be turned into a MATLAB script to calculate the largest LE of the fractional-order discrete-time SIR epidemic model with vaccination (11). Because it is difficult to predict the behavior of the system using analytical techniques, we must rely on approximate numerical approaches, which can be achieved by using MATLAB software. For further information on how to use numerical simulation of some chaotic systems using MATLAB, see the link in [53]. The results of the MLE for system with the same parameters and with initial conditions (S(0), I(0)) = (70, 30) are shown in Figure 7b. As can be seen, the system has positive LE, which shows that the MLE is compatible with the bifurcation diagram depicted in Figure 7a. We can also see that some values of Lyapunov exponents are negative, meaning that the system has periodic attractors, and this is illustrated further by the phase portraits plotted in Figure 5.

The 0−1 Test
To test chaotic behavior in nonlinear systems, one can use a 0−1 test method [54]. Unlike the Lyapunov exponent approach, the 0−1 test acts directly on the time sequence, therefore eliminating the necessity for phase reconstruction. In particular, this technique operates on finite points (S(n)) n=1,...,N and an arbitrary number c ∈ (0, π). According to the time series (S(n)), two terms referred to as the translation variables can be defined for n = 1, N as p c (n) = ∑ n k=1 S(k) cos(kc) and q c (n) = ∑ n k=1 S(k) sin(kc). We represent the mean square displacement by The asymptotic growth rate K c is given by the definition The asymptotic growth rate K = median(K c ) or the plotting of p c and q c on the p-q plane may be used to assess if chaos occurs on the fractional discrete SIR epidemic model. This means that the dynamics of the proposed fractional discrete SIR epidemic model are nonchaotic when K is close to 0 and the behavior of trajectory in the (p c -q c ) plane is bounded; on the other hand, when K is close to 1 and the trajectory in the (p c -q c ) plane exhibits Brownian-like behavior, the dynamics of the fractional discrete SIR epidemic model are chaotic.
The 0−1 test of the fractional discrete-time SIR epidemic model with vaccination (11) is applied directly to the state S(n). Figure 8 shows the (p c -q c ) plots of the fractional discrete SIR epidemic model (11) with commensurate fractional order. Observing the p c and q c trajectories, it is clear that the p c and q c trajectories display Brownian-like behavior when γ = 1 and γ = 0.94, indicating that the fractional SIR model is chaotic. Otherwise, when γ 3 = 0.99 and γ 3 = 0.966, the trajectories in the (p c -q c ) plane exhibit a bounded behavior, so the SIR model with commensurate order yields regular dynamics. Figure 8. The p-q trajectories of the 0−1 test of the fractional discrete SIR epidemic model (11) with commensurate orders (γ 1 = γ 2 ).

Incommensurate Fractional Order
In the following, the chaotic behavior of the proposed fractional SIR model with incommensurate fractional order (11) are carefully analyzed via the computation of bifurcation diagrams, Lyapunov exponents, phase portraits and 0−1 test method. Similarly, the effects of system parameter α and the fractional-order values on the dynamics of the model are illustrated in detail. We examine the case of incommensurate orders because it is more representative of reality than the case of commensurate orders, as each variable changes and moves independently of the others, implying that the rank of the influencer varies from one equation to another.

Bifurcation Diagram and Maximum LEs
Here, we discuss the dynamics of the fractional discrete SIR model (11) by varying α from 3 to 4.5 with the step size ∆α = 0.001. Figure 9 displays the bifurcation diagrams for the fractional order values (γ 1 , γ 2 ) = (0.94, 0.97), (γ 1 , γ 2 ) = (0.94, 0.99), and (γ 1 , γ 2 ) = (0.99, 0.94), respectively. From Figure 9, we see that the system exhibits periodic doubling or flip bifurcation. When α increases from 3, stability begins with a one-period orbit, then periodic orbits of periods 2, 4, and 8 are seen, which eventually evolve into chaos. In addition, we observe that the states of the system (11) are influenced by the fractionalorder values γ i and the system parameter α, as shown in Figure 9. As can be observed, when the incommensurate fractional orders (γ 1 , γ 2 ) are changed, the interval in which the chaos may be found shifts. Namely, when (γ 1 , γ 2 ) = (0.94, 0.97) the system is chaotic at α ∈]3.822, 3 (11) with incommensurate fractional orders shows a greater variety of chaotic attractors than the fractional discrete SIR epidemic model with commensurate fractional orders. To further investigate the effect of the incommensurate fractional orders on the dynamics behavior of system (11), Figure 10 displays the bifurcation diagram and the maximum Lyapunov exponents with γ 1 as a bifurcation parameter where the initial conditions are (S(0), I(0)) = (70, 30) , the system parameter is α = 4, and the fractional order is γ 2 = 0.94. We can see that the states of the system diverge towards infinity as soon as γ 1 goes below 0.8383 and when γ 1 increases, chaos occurs with certain periodic orbits. We can also see that when γ 1 approaches 1, the states become totally periodic. In Figure 10b, the maximum, Lyapunov exponents calculated for the same parameters and initial conditions as in Figure 10a are shown. Obviously, the maximum Lyapunov exponents have a positive number indicating that chaos exists, which agrees with the corresponding bifurcation diagram in Figure 10a. Figure 9. Bifurcation diagrams of the fractional discrete SIR model (11) for different incommensurate fractional-order values. Now, the dynamic behavior with the variation of the fractional order γ 2 is studied for γ 1 = 0.94. The bifurcation diagram and its corresponding maximum LE are illustrated in Figure 11. We can see that the dynamical behavior of the fractional discrete SIR epidemic model (11) evolves from periodic to chaos as γ 2 increases. In particular, the proposed SIR model is chaotic when γ 2 ∈]0.888, 0.9494[∪]0.951, 0.979[∪]0.9836, 1]; where the maximum LE is positive. The results illustrate that the dynamics behavior of the SIR model is affected by the order γ 2 . To illustrate the dynamics of the fractional discrete SIR epidemic model (11) better, phase portraits with different values of (γ 1 , γ 2 ) are presented in Figure 12. From Figure 12, we notice that the proposed SIR model shows different dynamic behaviors for these corresponding different fractional-order values. To get a better understanding of the influence of the fractional order on the SIR epidemic model, and in light of previous numerical findings, we compare the recent results obtained from the integer-order SIR model with the results obtained from the fractionalorder SIR model, which are presented in Table 1. It can be observed that the maximum number of susceptible cases obtained from the fractional-order SIR model is identical to the maximum number obtained from the integer-order SIR model, which is 70, whereas the minimum number of susceptible cases obtained from the fractional-order SIR model is less than the minimum number obtained from the integer-order SIR model. On the other hand, we see that the maximum and minimum numbers of infected cases predicted in the fractional-order SIR model are not identical to those expected in the integer-order SIR model. This confirms the effect of the fractional order on the SIR epidemic model in predicting the number of susceptible individuals and infected individuals. Similarly to the commensurate fractional orders, the 0−1 test was used to evaluate the fractional discrete SIR epidemic model (11) with incommensurate fractional order. The translation components p c and q c in the (p c -q c ) plan are illustrated in Figure 13. As can be observed, for (γ 1 , γ 2 ) = (0.91, 0.94) and (γ 1 , γ 2 ) = (0.94, 0.91), the trajectories p c and q c display Brownian-like behavior, and a bounded behavior for (γ 1 , γ 2 ) = (0.94, 0.983) and (γ 1 , γ 2 ) = (0.982, 0.94), which confirms that the fractional discrete SIR epidemic model (11) has chaotic attractors for (γ 1 , γ 2 ) = (0.97, 0.94) and (γ 1 , γ 2 ) = (0.97, 1), and has a regular behavior for (γ 1 , γ 2 ) = (0.99, 0.94) and (γ 1 , γ 2 ) = (0.97, 0.9). On the other side, Figure 14 depicts the asymptotic growth rate K with fractional orders (γ 1 , γ 2 ), in which α = 4 and (S(0), I(0)) = (70, 30). As we can see, for the majority of values of γ 1 and γ 2 , the asymptotic growth rate K approaches 1, implying that the fractional-order discrete SIR model with vaccination exhibits a chaotic behavior. These findings support very well the results of the bifurcation diagrams and the maximum Lyapunov exponents obtained before.  Figure 13. The p-q trajectories of the 0−1 test of the fractional discrete SIR epidemic model (11) with incommensurate orders (γ 1 , γ 2 ).

Complexity Analysis of the Fractional Discrete SIR Epidemic Model with Commensurate and Incommensurate Fractional Orders
For assessing the dynamic characteristics of chaotic systems, one technique is to consider the complexity of the chaotic characteristic. The model becomes more chaotic as the level of complexity increases. In this section, the approximate entropy and the C 0 complexity algorithm are used to assess the complexity of the fractional-order discrete-time SIR epidemic model with vaccination.

C 0 Complexity
Based on the inverse Fourier transform, the analysis of the complexity of chaotic systems was carried out using the C 0 algorithm. It divides the time series of the system into two components, a series of regular and a series of irregular parts, where the series of irregular parts is what we need. For a sequence [S(0), S(1), . . . , S(N − 1)] with a length of N, and a control parameter r, the algorithm process is defined as follows [55]. Firstly, the discrete Fourier transform of {S n } is determined by and the mean square value is calculated as The inverse Fourier transformation ofX N is defined as Then, the C 0 complexity is given by: The C 0 complexity of the fractional-order discrete-time SIR epidemic model with vaccination (11) with varying commensurate fractional-order value γ and incommensurate fractional-order values γ 1 , γ 2 are calculated and the result is shown in Figure 15. Interestingly, in the case of commensurate fractional-order values γ, as with the MLE, the C 0 complexity value of the fractional SIR model increases rapidly when γ i decreases. On the other hand, the fractional discrete SIR epidemic model (11) with incommensurate fractionalorder values, in contrast to the case when the model has commensurate fractional-order values, has more complexity when γ 2 approaches 1. Thus, we can see that the C 0 algorithm can measure the complexity effectively. Figure 15a illustrates that model (11) has a higher complexity when γ ∈ (0.9712, 0.9082]. When γ 1 = 0.94, the high complexity region in Figure 15c also exists in the range of γ 2 ∈ [0.9836, 1], and the C 0 complexity increases with the increase of fractional order γ 2 .

Approximate Entropy (ApEn)
The approximate entropy is a measure of regularity that quantifies the level of complexity within systems generated by a time series. Generally, time series with larger values of ApEn are considered to be more complex [56]. Note that, the approximate entropy value is dependent on two essential parameters: the similarity tolerance r and the embedding dimension m. In this article, we took m = 2 and r = 0.2std(S) where std(S) represents the standard deviation of the data S. Theoretically, the ApEn is calculated as follows: where Φ m (r) is denoted by where C m i (r) = K n−m+1 , i ∈ [1, n − m + 1], K is the number of X(j) such that d(X(j), X(i)) ≤ r, and X(i) = [x(i), x(i + 1), ..., x(i + m − 1)].
The ApEn complexity of the fractional-order discrete-time SIR epidemic model with vaccination (11), with both commensurate and incommensurate orders (γ 1 , γ 2 ) was ana-lyzed and the findings are reported in Table 2. It can be seen that the complexity of the fractional discrete SIR epidemic model (11) varies when γ i (i = 1, 2) varies, and the highest ApEn is found when the model is chaotic, which agrees very well with the maximum LE results. As a result, in order to obtain a relatively high structural complexity, we must be cautious when choosing the values of γ i in system (11).

Conclusions
In this paper, we dealt with the dynamics of a new fractional-order discrete-time SIR epidemic model with vaccination with commensurate and incommensurate orders. Through phase portraits, bifurcation diagrams, and maximum LEs, the complex dynamics of the proposed system were discussed. In addition, we also calculated the 0−1 test, approximate entropy (ApEn), and C 0 complexity of the fractional-order discrete SIR epidemic model for different fractional order values, all of which intended to demonstrate and quantify the complex dynamics of the system. Results indicated that the fractionalorder model is more complex than the integer-order model. Furthermore, we showed that the incommensurate fractional orders have a greater effect on the behavior of the system than when the fractional orders are commensurate. The reasonable range of commensurate fractional orders is between γ = 0.8712 and γ = 1 while the reasonable range of incommensurate fractional orders is between γ 2 = 0.77 and γ 2 = 1. Throughout this work, numerical simulations were used to explain all the findings. Due to the current spread of the COVID-19 epidemic and the onset of the implementation of various vaccination strategies in various countries, we will attempt to use the fractal properties of most of the COVID-19 series in order to establish a connection between the proposed fractional model and the stochastic self-affine characteristics of the COVID-19 times series in future research. In addition, the determination of the coupled time-fractional differential equations of the proposed fractional model using the Caputo fractional derivative, as well as the study of some complex epidemiological models will be considered in our future research.