On the Dynamical Complexity of a Seasonally Forced Discrete SIR Epidemic Model with a Constant Vaccination Strategy

1 School of Mathematics, Iran University of Science and Technology, Narmak, Tehran, Iran 2Instituto Superior de Engenharia de Lisboa (ISEL), Department of Mathematics, Rua Conselheiro Emı́dio Navarro 1, 1949-014 Lisboa, Portugal 3Center for Mathematical Analysis, Geometry and Dynamical Systems, Mathematics Department, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal 4Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal


Introduction
Understanding the dynamical behaviors of emergent infectious diseases in humans is viewed with increasing importance in epidemiology.In particular, mathematical modeling has been used to provide useful insights to enhance our comprehension of complex processes associated with the pathogenesis of diseases, as well as to quantify the likely effects of different intervention/control strategies (see, for example, [1,2]).In fact, epidemiological models are now widely used as more epidemiologists realize the role that modeling can play in basic understanding of complex factors and policy development.More precisely, as far as infectious diseases are concerned, very recent works have been conducted in the literature characterizing different transmission processes and providing new insights for elimination/control, since these diseases impose a considerable burden on human population ( [3][4][5]).In addition, useful discussions have taken place on the disease spread with space.In many cases the spatial structure of the population cannot be ignored.Consequently, the study of the spread of diseases in both time and space revealed different types of spatial patterns.Reviews on the transitions between patterns and their underlying mechanisms have been proposed in [6,7].
Generally, in the theory of epidemics, there are two kinds of mathematical models: the continuous-time models described by differential equations and the discrete-time models described by difference equations.Mathematical 2 Complexity models (continuous and discrete) describing the population dynamics of infectious diseases have been playing an important role in enhancing our understanding of complex epidemiological patterns driven by seasonality.The continuoustime epidemic models described by differential equations have been widely investigated in the literature.In recent years, due to infectious diseases data collected based on periods of days, weeks, months, and years, we have seen more attention being given to a great deal of practical real life problems in epidemiology depicted by discrete-time epidemic models described using difference equations.Although, in recent times, an increasing attention has been given to the discrete epidemic models (see [2,8] and the references cited therein), the research in the context of the nonlinear dynamic characteristics of the discrete systems is relatively few, compared with the amount of studies performed in the literature about the continuous systems.With their unique dynamic characteristics, the discrete systems are able to depict many practical problems in the real life, in particular in the context of epidemiology (see [9][10][11][12] and the references cited therein) The advantages of the discrete-time approach are multiple in the study of epidemic models ( [13,14]).Firstly, with statistical epidemic data compiled from given time intervals and not continuously, the difference models appear to be more realistic and accurate to describe epidemics.A second reason is the fact that the discrete-time models can be used as natural simulators of the continuous cases.As a consequence, we can not only analyze with remarkable accuracy the dynamical behavior of the continuous-time models, but also assess the effect of larger time steps.At last, the use of discrete-time models makes it possible to use a combination of methods from the theory of nonlinear dynamical systems for the study of mappings (and lattice equations), from the integrability and/or chaos points of view.The particularly rich dynamical behaviors of discrete-time models, involving bifurcations and chaos, is one of the noteworthy and eye-catching features of their dynamics.
As pointed out in paper [15], there are mainly two ways to construct discrete-time epidemic models: (i) by making use of the property of the epidemic disease (see [16,17]) and (ii) by discretizing a continuous-time epidemic model using techniques such as the forward Euler scheme and the Mickens nonstandard discretization (see [15]).Let us consider the following classical continuous-time Susceptible-Infected-Recovered (SIR) forced epidemic model described by the differential equations The susceptible group refers to individuals with density () who have never come into contact with the disease at time , but they could catch it.Infectious individuals are capable of spreading the disease to those in the susceptible category and remain in the infectious compartment until their recovery, with density ().The recovered class refers to individuals with density () who have been infected and then recovered from the disease.The recovered individuals are immune for life and are not able to transmit the infection to others.They are essentially removed from the population and play no further role in the dynamics.Epidemics are continuously fueled by the constant supply of new susceptibles that arise due to the birth of new individuals.
The vital parameters of this system are , the birth rate, and , the death rate.The parameter  denotes the recovery rate.Under seasonality, a commonly used scheme for the contact rate (or transmission rate) () takes where  0 gives the mean contact rate (the base of transmission rate).All the parameters of system (1) that have been described, i.e., , , , and  0 , are nonnegative and smaller than 1.The parameter , with 0 ≤  ≤ 1, represents the strength of the seasonal forcing (measuring the degree of seasonality);  is a -periodic function of zero mean and  is scaled in units of years.In certain studies,  is considered as a sinusoidal function [18].Nevertheless, the transmission rate of many childhood infections and seasonal influenza abruptly changes by the host behavior (cycles of human aggregation such as migrations and school holidays).Inspired by the approach adopted in [19], we assume that  can be expressed as In this context, we have () =  0 (1 + ) in the high season and () =  0 (1 − ) in the low season.For simplicity and clarity reasons, we have assumed that the lengths of the high and low seasons are the same.The newborns are susceptible individuals and an exposure of a susceptible to an infective is an encounter in which the infection is transmitted.In this context, the contact rate, (), is the average density of susceptible in a given population contacted per infective individuals per unit of time.Therefore, the product ()() denotes the rate of the total density of susceptible infected by one infective and ()()() represents the rate of infection of the susceptible by all infectives.
Applying the forward Euler scheme to model (1), we obtain the following discrete-time forced SIR epidemic model where ℎ is the step size.In this paper, we focus on the dynamical behaviors of model (4), which include chaos phenomenon caused by the change of the time step ℎ.The first two equations of the model (4) do not include   and the third equation is a linear equation of   .Thus, we are going to focus on variables (  ,   ) and the dynamical behaviors of model which only includes   and   .
Seasonal forces in epidemic systems shape the population dynamics, particularly the spread of the infectious diseases ( [20,21]).We will see that the variation of the time step ℎ can result in deterministic chaos, even considering a constant rate, () =  0 .In the following sections of our study, the time step ℎ is chosen as a bifurcation parameter to study different complex dynamical behaviors of model (5) in two regimes: epidemics without seasonality ( = 0), taking () =  0 , (Section 2) and epidemics with seasonal contact rate ( > 0), taking () =  0 (1 + ()), with () given by (3) (Section 3).
In our analysis, if not otherwise specified, we will set  = 0.72,  = 0.4,  = 0.24,  0 = 0.445.A variety of interesting dynamical behaviors will occur for values of the strength of the seasonal forcing, , such that  ∈ [0, 0.0075] and for the time step ℎ as a control parameter, with ℎ ∈ [2.92, 3.4], also an interval for which the rich dynamics takes place.In fact, it is important to emphasize that the time step size ℎ is taken as a bifurcation parameter in order to detect different complex dynamics.In particular, by choosing ℎ as a bifurcation parameter in the considered interval, [2.92, 3.4], we are able to exhibit a variety of dynamical behaviors from period 1 and period-doubling bifurcations, passing through interesting rich qualitative dynamics (flip bifurcation, Hopf bifurcation) to chaos.This way, the comprehensive range of [2.92, 3.4] for values of ℎ was chosen to fulfil the aim of displaying a window of global noticeable qualitative features of the dynamics of the studied discrete-time epidemic model.In all of our numerical simulations, we take as initial values of the dynamical variables  0 = 1.2 and  0 = 0.4, displaying, in a variety of situations, their asymptotic dynamical behaviors.

Epidemics without Seasonality
Taking the discrete-time forced epidemic model ( 5), we consider in this section a constant contact rate (or constant transmission rate).In this context, after the identification of global noteworthy features of the dynamics and the study of the maximal Lyapunov exponent, a biologically meaningful family of unimodal-type maps, associated with the density of susceptible, is detected varying the time step ℎ.Based on the theory of symbolic dynamics, we are going to compute an important numerical invariant associated with these iterated maps, the topological entropy.

Global Noticeable Qualitative Features of the Dynamics:
The Maximal Lyapunov Exponent.In order to start gaining some insights into the qualitative dynamic behavior of model ( 5), we display in Figure 1(a) the (ℎ,   ,   )-bifurcation diagram and in Figure 1(b) the (ℎ,   )-bifurcation diagram as a result of the variation of the time step ℎ, giving a particular emphasis to the dynamical variable   .When the time step ℎ is sufficiently small, the dynamical behaviors of the discretetime model ( 5) are similar to the ones exhibited by the continuous-time SIR epidemic model.At increasing values of ℎ, complex behaviors appear, namely, a period-doubling route to chaos (the variable   has a similar dominant qualitative behavior to   ).The threshold value of the time step ℎ that causes the bifurcations and chaos can be effectively computed according to the method developed in [10].Despite its importance, the explicit computation of this threshold value for ℎ goes beyond the central objectives of the present study.
In the context of epidemic models, it is usual to define the so-called basic reproduction number,   .In the case of constant contact rate, () =  0 , it is given by which denotes the average number of secondary infections generated by an initial population of infected individuals over their lifetimes.This number can be interpreted as a "disease-spreading impact" of one infected individual in the population with only susceptible individuals.
From the literature on the existence of nonnegative equilibria of discrete-time SIR type models (see, for instance, [10,11]), the following are known: (i) Model ( 5) always has a disease-free equilibrium (DFE) corresponding to the population with no infected individuals.In this case, the epidemic cannot maintain itself in the population and the DFE point ( * 0 ,  * 0 ) is locally stable.(ii) If   > 1, then model (5) also has an endemic equilibrium (EE) corresponding to a self-sustaining group of infected individuals in the population.In this case, the EE point ( * 1 ,  * 1 ) is locally stable.As a consequence, the disease will remain permanently endemic in the population.
A powerful numerical tool to examine the degree of chaoticity of the discrete-time model ( 5) is the computation, and study of the variation, of the maximal Lyapunov exponent, as a function of the time step ℎ.The largest Lyapunov exponent is the average growth rate of an infinitesimal state perturbation along a typical orbit.In particular, this numerical invariant is a convenient indicator of the exponential divergence of initially close points, characteristic of the chaotic attractors (please see [22][23][24]).A discussion about the Lyapunov exponents as a quantitative measure of the rate of separation of infinitesimally close trajectories, as well as a computation method, can be found in [25].
Let us consider the following nonlinear map: where F : R  → R  is a differentiable vector function and y is a  vector.Let DF(y) denote the matrix (  /  ) of partial derivatives of the components   at y. Then the corresponding matrix of partial derivatives for the  ℎ iterate F  of F is given by Considering we have Then In this context,   can be decomposed into a product of an orthogonal matrix and an upper triangular matrix with positive diagonal elements.More precisely, where   ,  +1 ,   , and  +1 are the upper triangular matrices with positive diagonal elements.Taking { ()  ,  = 1, . . ., } as the set of diagonal elements of   , the Lyapunov exponents are given by with  = 1, . . ., .
We display in Figure 1(c) the variation of the maximal Lyapunov exponent of system (5),  max , with ℎ along the interval [2.92, 3.4].In agreement with the bifurcation diagrams, we have  max < 0 within the periodic windows and  max = 0 at the bifurcation points.The positivity of the maximal Lyapunov exponent is considered one of the criteria of chaotic motions.For values of ℎ in the interval [3.014, 3.182], a phenomenon of homeochaos was found corresponding to 0 <  max ≪ 1.
The region in red displayed in Figure 1(b) and in Figure 1(c), corresponding to values of the time step ℎ in the interval  = [  ,   ] = [3.25500,3.31097], is associated with salient, eye-catching, and noteworthy features of the dynamics of the variable   that will be studied in the following section.

Iterated Maps, Symbolic Dynamics, and Topological
Entropy.In this paragraph, a considerable relevance is given to the dynamical variable   studied for ℎ within the interval  = [  ,   ] = [3.25500,3.31097], presented previously.We show in Figure 2(a) a typical -phase diagram and in Figure 2(b) the corresponding variation of the iterates of   for ℎ = 3.31075.In order to gain additional significant qualitative insights into the principles and mechanisms underlying the dynamics of these iterates of   , we started our study by recording the successive relative (local) maxima of   .
Each local maximum corresponds to a peak in the susceptible series of iterates.Considering the pairs ( max  ,  max +1 ), where  max  denotes the  ℎ local maximum, an unimodaltype iterated map emerges.As shown in Figure 2(c) the data from the series of values of   , corresponding to successive points ( max  ,  max +1 ), appear to fall on a logistic curve.Indeed, treating the graph as a representation of a function  max +1 = ( max  ) allows us to reveal particularly interesting features of the dynamics.The obtained iterated maps dynamically behave like unimodal maps, which have found significant applications on symbolic dynamics theory, mapping an interval  = [, ] into itself.
In the next lines, we will explain, following closely the work carried out in [26] in the context of symbolic dynamics, the procedure used to compute an accurate estimate of a numerical invariant that arises as a measure of the amount of chaos, the topological entropy.This numerical invariant is a central measure related to orbit growth ( [27,28]).The field of symbolic dynamics evolved as a tool for analyzing general dynamical systems by taking the state space and considering its partition into a finite number of regions, each of which labeled with a given symbol.We obtain a symbolic trajectory Complexity by writing down the sequence of symbols corresponding to the successive partition elements visited by a point following some trajectory in the phase space.The symbolic coding of the intervals of a piecewise monotonic map and the study of these symbolic sequences allow us to analyze qualitative aspects of the dynamical system in two perspectives: the kneading theory and the theory of Markov partitions.In this work, we follow the results concerning Markov partitions ( [29][30][31]).
Let us consider a unimodal map  and the interval  subdivided into the sets   = [, [,   = {}, and   =], ].The function  is monotonically increasing for  ∈   and monotonically decreasing for  ∈   (Figure 2(c)).We call a lap of  each of such maximal intervals where the map  is strictly increasing or strictly decreasing.The total number of distinct laps is called the lap number of  and it is usually denoted by ℓ = ℓ().The left and the right subintervals are labeled by the letters  and , respectively, and the set   will be denoted by .The symbolic sequence starting from () plays an important role in the symbolic dynamics of one-dimensional maps and it is called kneading sequence.Consequently, let us consider the orbit of the critical point of , (), obtained by iterating the map, that is, () = {  :   =   (),  ∈ N}.From this numerical orbit, (), we get a symbolic sequence (symbolic orbit)  =  1  2 . . .  . .., where that characterizes the dynamics.When () is a -periodic orbit, we obtain a sequence of symbols that can be characterized by a block of length , the kneading sequence  () =  1  2 . . . −1 .The orbit (), which is made of ordered points   , determines a partition P (−1) of the invariant range  = [ 2 (), ()] = [ 2 ,  1 ] into a finite number of subintervals labeled by  1 ,  2 , . . .,  −1 .This partition is associated with a ( − 1) × ( − 1) transition matrix  = [  ], where the rows and columns are labeled by the subscript of subintervals and the matrix elements are defined as Now, taking up the generalized problem of characterizing the chaoticity of the dynamics with symbolic dynamics theory, the topological entropy represents the exponential growth rate for the number of orbit segments distinguishable with arbitrarily fine but finite precision.The topological entropy describes in a suggestive way the exponential complexity of the orbit structure with a single nonnegative real number ( [28]).More precisely, the growth rate of the lap number of   ( ℎ -iterate of ) is and the topological entropy of a unimodal interval map , denoted by   (), is given by where  max (()) is the spectral radius of the transition matrix () ( [29,30,32]).The following example illustrates the computation of the topological entropy using the previously established procedure.
Example.Let us consider the orbit of a turning point defined by the period-6 kneading sequence () ∞ .Putting the orbital points in order we obtain The corresponding transition matrix is which has the characteristic polynomial The growth number () (the spectral radius of matrix ()) is 1.96595 . ... Therefore, the value of the topological entropy can be given by The study of the kneading sequences allows us to identify symbolic periodic orbits.In particular, with periods  ≤ 5, we have the following: and 1-period - ∞ .The identified kneading sequences correspond to logistic-type iterated maps with different levels of complexity reflected in different values of the topological entropy.Table 1 gives us information about some nuclear kneading sequences studied in terms of its characteristic polynomial and its topological entropy.
The variation of the topological entropy in the parameter region depicted in Figure 2(e) is in agreement with the bifurcation diagram of Figure 2(d).As we can observe, the dynamics of the system associated with positive topological entropy highlights a region of the parameter space where chaos occurs.Our results suggest that the dynamics is very sensitive to the time step ℎ.The variation of higher values of this parameter has a profound and marked effect on the dynamical behavior of the system, where a transition from order to chaos takes place.The positive topological entropy, which increases for growing ℎ, allows us to identify the chaotic regimes.As a consequence, the feature of the original model that we are studying-the temporal dynamical behavior of the successive local maxima of the susceptible individuals,   -is associated with regimes of chaotic behavior.
It is fundamental to emphasize that the existence of onedimensional unimodal-type iterated maps is a particularly remarkable feature of the dynamics without seasonality (i.e., with constant contact rate).With the identification of these maps, our attention is directed to a striking dynamical characteristic of the model without seasonality-a local maximum  max  of the densities of susceptible predicts the following pick  max +1 of the same population.
As pointed out before, the data from the chaotic time series (Figure 2(b)) appears to fall neatly on a curve (Figure 2(c)) (notice that there is almost no thickness to the graph).By this trick, we were able to extract order from chaos.From the biological point of view, treating the graph of Figure 2(c) as a plausible representation of a function  max +1 = ( max  ) allows us to tell a lot about the dynamics: given  max 0 , we can predict max 1 by  max 1 = ( max 0 ), then use that information to predict  max 2 = ( max 1 ), and so on, bootstrapping our way forward in time by iteration.
In the following sections, we are going to point out interesting characteristics of the dynamics that arise with the introduction of a degree (or amplitude) of seasonality .

Epidemics with Seasonal Contact Rate
For the epidemics with a seasonal contact rate, the different chaotic regimes are mainly studied through bifurcation diagrams and the computation of the maximal Lyapunov exponent.This analysis includes the construction of a 3Dplot that illustrates the variation of the maximal Lyapunov exponent within the (, ℎ) parameter space.
To start with, we exhibit in Figure 3 the graph of the step function () that we are going to use in the model.With the same aim of retaining noticeable features of the dynamics under seasonality, we take now  = 0.002 and we show in Figure 4(a) the (ℎ,   ,   ) bifurcation diagram and in Figure 4(b) the (ℎ,   ) bifurcation diagram, as a consequence of the variation of the time step ℎ.With  > 0, a nonautonomous system emerges.The differences between Figures 4 and 1, of the previous section, are a signature of the presence of a term in the system directly depending on time .As suggested by Figure 4, the introduction of the degree (or amplitude) of seasonality, , is responsible for inducing complexity into the population dynamics.In Figure 4(c), we show the corresponding variation of the maximal Lyapunov exponent,  max , with the time step ℎ, giving rise to a particular set of maximal Lyapunov exponents, denoted by .With the 3D graph of Figure 5, we intend to illustrate the variation of the maximal Lyapunov exponent in the (ℎ, ) parameter space.In this context, for each  we have a particular set , obtained with the variation of the time step ℎ.As expected in this regime with seasonal fluctuations, a region for which  max > 0 (a stated signature of chaotic behavior) is obtained.The variation of the maximal Lyapunov exponent also allows us to identify a window where a phenomenon of homeochaos occurs, corresponding in this case to ℎ ∈ [3.14, 3.17] and characterized by 0 <  max ≪ 1.
It is known that the onset of chaos can cause the population to run a higher risk of extinction due to the unpredictability of the dynamics.In this context, the density of the infected may be out of control.However, in the real world, the density of infected needs to be under control; otherwise it will be harmful to the health of people worldwide.As a consequence, how to control chaos in the epidemic model is a very important aspect to study, which motivates the investigation conducted in the following section.We will see that a planned constant vaccination acts as an effective control strategy.

Dynamics of the Model under a Constant Vaccination Strategy
In this section, we carry out a study of the discrete SIR epidemic model with seasonal fluctuations under a planned vaccination regime.According to the conventional constant vaccination strategy, all newborn infants should be vaccinated, and V 0 is the proportion of those vaccinated successfully (with 0 < V 0 < 1).Indeed, constant vaccination takes action in the first equation of model (5), reducing the birth rate, , of susceptible, so that the system becomes where the term (1−V 0 ) = −V 0  incorporates the density of newborn infants vaccinated, V 0 .In order to gain some insights into the possible effect of the constant vaccination regime in the dynamics, we start our analysis by considering the mean value of the positive maximal Lyapunov exponents.With this purpose, let us denote a positive maximal Lyapunov exponent by  + max ∈  + , with  + the set of all values  + max that result from the variation of the time step ℎ, for fixed values of  and V 0 .In Figure 6(a) we represent the behavior of the mean value of the positive Lyapunov exponents,  + max , as a function of the proportion of susceptibles who received constant vaccination, V 0 , for a lower degree of seasonality and also for a higher degree of seasonality .The respective derivatives of  + max in order to V 0 ,  + max /V 0 , are depicted in Figure 6(b).A prominent feature emerges from these two graphs of Figure 6: the value of V 0 = 0.02 seems to act as a threshold, having a decisive effect on the dynamics.In this context, for values V 0 < 0.02 the dynamics is still exhibiting chaotic behavior, with  max > 0 for some time steps corresponding to the points inside the dashed circle of Figure 7 As a consequence, with the chaotic behavior suppressed, the system converges to controlled population densities.This threshold property of V 0 is illustrated in Figure 7 taking, as an example,  = 0.002.All our numerical simulations for other values of  also confirm V 0 = 0.02 as a threshold of the dynamics (results not shown).

Final Considerations
In this work, we have converted the classical SIR epidemic model with seasonal fluctuations into a discretized system and discussed its dynamical behaviors.The discrete model can result in a much richer set of patterns than the corresponding continuous-time model and it is more effective in practice.Taking the time step ℎ as a bifurcation parameter, the discrete model displays particularly interesting dynamical behaviors, such as period-7 orbits, period-doubling cascades, and chaotic sets, where susceptible and infective coexist.
Having started with the analysis of the dynamical features of the epidemics without seasonality, a family of unimodaltype iterated maps has been identified, considering the evolution of the local maxima of the susceptible density.The existence of these applications offered us an appropriate dynamical context to use symbolic dynamics to compute the topological entropy.The positivity of this important numerical invariant revealed regimes of chaotic behavior associated with the temporal dynamical evolution of the successive local maxima of the susceptible individuals,   .Noteworthy for its biological meaning, the one-dimensionality of the identified iterated maps implies that we can predict, with a remarkable accuracy, the next pick of the population density of susceptible based on its previously recorded local maximum.
In the context of epidemiology, more realistic dynamics may be achieved by taking into account sources of seasonal variation in the contact rate attributed to social behavior, namely, the schedule of school year and seasonal changes in weather conditions.The introduction of a seasonal forcing provides the discrete SIR model with the potential to generate more complex oscillatory and chaotic dynamical regimes.Given the importance of controlling this chaotic behavior, i.e., the necessity of having predictable densities for the epidemic populations, we have applied a planned constant vaccination of the newborn infants, which has acted as an effective control strategy.More precisely, the chaotic epidemic outbreaks, which appeared as a result of the seasonal variations in the contact rate, have been suppressed by the used vaccination scheme.In particular, we have identified a specific threshold value, for the proportion of newborn infants vaccinated successfully, which seems to be independent of the degree of the seasonality.This study illustrates how an approach integrating theoretical reasonings and numerical experiments can contribute to our understanding of important biological models and provide trustworthy explanation of complex phenomena witnessed in biological systems.

Figure 1 :
Figure 1: (a) 3D bifurcation diagram corresponding to the dynamical variables   and   , with the time step ℎ as the control parameter.(b) Dynamics of the variable   , taking also ℎ as the bifurcation parameter.(c) Variation of the maximal Lyapunov exponent of the system (5) with ℎ.In all the previous situations () =  0 = 0.445 and ℎ ∈ [2.92, 3.4].The displayed region corresponds to the interval of values ℎ, ℎ ∈  = [  ,   ] = [3.25500,3.31097], and is associated with noticeable features of the dynamical behavior of   , as we will see in our study below.

Figure 4 :Figure 5 :
Figure 4: (a) 3D bifurcation diagram corresponding to the dynamical variables   and   , with the time step ℎ as the control parameter.(b) Dynamics of the variable   , taking also ℎ as the bifurcation parameter.(c) Variation of the maximal Lyapunov exponent of the system (5) with ℎ.In all the previous situations  = 0.002 and ℎ ∈ [2.92, 3.4].
(a).For values V 0 ≥ 0.02 those types of points tend to disappear (please see the circles of Figures7(b) and 7(c)).