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

Authors: Jalil Rashidinia; Mehri Sajjadian; Jorge Duarte; Cristina Januário; Nuno Martins
Journal: Complexity
 (2018)

---

## Abstract

In this article, we consider the discretized classical Susceptible-Infected-Recovered (SIR) forced epidemic model to investigate the consequences of the introduction of different transmission rates and the effect of a constant vaccination strategy, providing new numerical and topological insights into the complex dynamics of recurrent diseases. Starting with a constant contact (or transmission) rate, the computation of the spectrum of Lyapunov exponents allows us to identify different chaotic regimes. Studying the evolution of the dynamical variables, a family of unimodal-type iterated maps with a striking biological meaning is detected among those dynamical regimes of the densities of the susceptibles. Using the theory of symbolic dynamics, these iterated maps are characterized based on the computation of an important numerical invariant, the topological entropy. The introduction of a degree (or amplitude) of seasonality,ε, is responsible for inducing complexity into the population dynamics. The resulting dynamical behaviors are studied using some of the previous tools for particular values of the strength of the seasonality forcing, ε. Finally, we carry out a study of the discrete SIR epidemic model under a planned constant vaccination strategy. We examine what effect this vaccination regime will have on the periodic and chaotic dynamics originated by seasonally forced epidemics.


---

## Body

## 1. 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–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 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 continuous-time 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–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 theforward Euler scheme and theMickens nonstandard discretization (see [15]). Let us consider the following classical continuous-time Susceptible-Infected-Recovered (SIR) forced epidemic model described by the differential equations(1)dSdt=α-μS-βtSIdIdt=βtSI-γ+μIdRdt=γI-μR.Thesusceptible group refers to individuals with density S(t) who have never come into contact with the disease at time t, 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 I(t). Therecovered class refers to individuals with density R(t) 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α, thebirth rate, and μ, thedeath rate. The parameter γ denotes therecovery rate. Under seasonality, a commonly used scheme for thecontact rate (ortransmission rate) β(t) takes(2)βt=β01+εϕt,where β0 gives themean 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 thestrength of the seasonal forcing (measuring the degree of seasonality); ϕ is a T-periodic function of zero mean and t 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(3)ϕt=-1if0≤t<T11ifT1≤t<2T1≕T,whereT1>0.In this context, we have βt=β01+ε in the high season and βt=β01-ε 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,β(t), is the average density of susceptible in a given population contacted per infective individuals per unit of time. Therefore, the product β(t)S(t) denotes the rate of the total density of susceptible infected by one infective and β(t)S(t)I(t) represents the rate of infection of the susceptible by all infectives.Applying theforward Euler scheme to model (1), we obtain the following discrete-time forced SIR epidemic model(4)Sn+1=Sn+hα-μSn-βnhSnInIn+1=In+hβnhSnIn-γ+μInRn+1=Rn+hγIn-μRn,where h 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 h. The first two equations of the model (4) do not include Rn and the third equation is a linear equation of Rn. Thus, we are going to focus on variables Sn,In and the dynamical behaviors of model(5)Sn+1=Sn+hα-μSn-βnhSnInIn+1=In+hβnhSnIn-γ+μIn,which only includes Sn and In.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 h can result in deterministic chaos, even considering a constant rate, β(t)=β0. In the following sections of our study, the time step h is chosen as a bifurcation parameter to study different complex dynamical behaviors of model (5) in two regimes:epidemics without seasonality(ε=0), taking β(t)=β0, (Section 2) andepidemics with seasonal contact rate(ε>0), taking β(t)=β01+εϕt, with ϕt 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 h as a control parameter, with h∈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 h is taken as a bifurcation parameter in order to detect different complex dynamics. In particular, by choosing h 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 h 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 S0=1.2 and I0=0.4, displaying, in a variety of situations, their asymptotic dynamical behaviors.


## 2. 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 h. Based on the theory of symbolic dynamics, we are going to compute an important numerical invariant associated with these iterated maps, the topological entropy.


### 2.1. 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 h,Sn,In-bifurcation diagram and in Figure 1(b) the h,Sn-bifurcation diagram as a result of the variation of the time step h, giving a particular emphasis to the dynamical variable Sn. When the time step h is sufficiently small, the dynamical behaviors of the discrete-time model (5) are similar to the ones exhibited by the continuous-time SIR epidemic model. At increasing values of h, complex behaviors appear, namely, a period-doubling route to chaos (the variable In has a similar dominant qualitative behavior to Sn). The threshold value of the time step h 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 h goes beyond the central objectives of the present study.Figure 1
(a) 3D bifurcation diagram corresponding to the dynamical variablesSn and In, with the time step h as the control parameter. (b) Dynamics of the variable Sn, taking also h as the bifurcation parameter. (c) Variation of the maximal Lyapunov exponent of the system (5) with h. In all the previous situations β(t)=β0=0.445 and h∈2.92,3.4. The displayed region corresponds to the interval of values h, h∈M=mi,mf=3.25500,3.31097, and is associated with noticeable features of the dynamical behavior of Sn, as we will see in our study below.
(a) (b) (c)In the context of epidemic models, it is usual to define the so-calledbasic reproduction number, Rp. In the case of constant contact rate, β(t)=β0, it is given by (6)Rp=β0μ+γ,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) (7)S0∗,I0∗=αμ,0,corresponding to the population with no infected individuals. In this case, the epidemic cannot maintain itself in the population and the DFE point S0∗,I0∗ is locally stable.(ii) IfRp>1, then model (5) also has an endemic equilibrium (EE) (8)S1∗,I1∗=γ+μβ0,-μβ0+αγ+μ,corresponding to a self-sustaining group of infected individuals in the population. In this case, the EE point S1∗,I1∗ 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 h. 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–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:(9)yn+1=Fyn,where F:Rm→Rm is a differentiable vector function and y is a m vector. Let DFy denote the matrix ∂Fi/∂yi of partial derivatives of the components Fi at y. Then the corresponding matrix of partial derivatives for the nth iterate Fn of F is given by (10)∂Fni∂yj=DFFn-1y⋯DFFyDFy.Considering(11)An=DFFn-1y⋯DFFyDFy,we have(12)An=DFyn-1⋯DFy1DF0.Then(13)An+1=DFynAn.In this context,An can be decomposed into a product of an orthogonal matrix and an upper triangular matrix with positive diagonal elements. More precisely,(14)Qn+1Rn+1=DFynQnRn,where Qn, Qn+1, Rn, and Rn+1 are the upper triangular matrices with positive diagonal elements. Taking Rijn,i=1,…,m as the set of diagonal elements of Rn, the Lyapunov exponents are given by (15)λi=limn→∞ln⁡Rijnn,with i=1,…,m.We display in Figure1(c) the variation of the maximal Lyapunov exponent of system (5), λmax, with h 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 h in the interval 3.014,3.182, a phenomenon of homeochaos was found corresponding to 0<λmax≪1.The region in red displayed in Figure1(b) and in Figure 1(c), corresponding to values of the time step h in the interval M=mi,mf=3.25500,3.31097, is associated with salient, eye-catching, and noteworthy features of the dynamics of the variable Sn that will be studied in the following section.


### 2.2. Iterated Maps, Symbolic Dynamics, and Topological Entropy

In this paragraph, a considerable relevance is given to the dynamical variableSn studied for h within the interval M=mi,mf=3.25500,3.31097, presented previously. We show in Figure 2(a) a typical SI-phase diagram and in Figure 2(b) the corresponding variation of the iterates of Sn for h=3.31075. In order to gain additional significant qualitative insights into the principles and mechanisms underlying the dynamics of these iterates of Sn, we started our study by recording the successive relative (local) maxima of Sn.Figure 2
(a)SnIn-phase diagram (h=3.31075). (b) Series of the iterations of Sn. (c) Unimodal-type iterated map corresponding to the successive points Snmax,Sn+1max, where Snmax denotes the nth local maximum (h=3.31075). (d) Bifurcation diagram of Snmax, taking h∈M=3.25500,3.31097 as the control parameter. € Variation of the topological entropy corresponding to the dynamics of Snmax for h∈M=3.25500,3.31097.
(a) (b) (c) (d) (e)Each local maximum corresponds to a peak in the susceptible series of iterates. Considering the pairsSnmax,Sn+1max, where Snmax denotes the nth local maximum, an unimodal-type iterated map emerges. As shown in Figure 2(c) the data from the series of values of Sn, corresponding to successive points Snmax,Sn+1max, appear to fall on a logistic curve. Indeed, treating the graph as a representation of a function Sn+1max=fSnmax 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 I=a,b 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 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–31]).Let us consider a unimodal mapf and the interval I subdivided into the sets IL=[a,c[, IC={c}, and IR=]c,b]. The function f is monotonically increasing for x∈IL and monotonically decreasing for x∈IR (Figure 2(c)). We call a lap of f each of such maximal intervals where the map f is strictly increasing or strictly decreasing. The total number of distinct laps is called the lap number of f and it is usually denoted by l=lf. The left and the right subintervals are labeled by the letters L and R, respectively, and the set IC will be denoted by C. The symbolic sequence starting from f(c) 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 f, Oc, obtained by iterating the map, that is, Oc=xi:xi=fic,i∈N. From this numerical orbit, Oc, we get a symbolic sequence (symbolic orbit) S=S1S2…Sj…, where(16)Sj=Liffjx<c,Sj=Ciffjx=cSj=Riffjx>c,that characterizes the dynamics. When O(c) is a k-periodic orbit, we obtain a sequence of symbols that can be characterized by a block of length k, the kneading sequence S(k)=S1S2…Sk-1C. The orbit Oc, which is made of ordered points xi, determines a partition P(k-1) of the invariant range I=f2(c),f(c)=x2,x1 into a finite number of subintervals labeled by I1,I2,…,Ik-1. This partition is associated with a (k-1)×(k-1) transition matrix M=[aij], where the rows and columns are labeled by the subscript of subintervals and the matrix elements are defined as aij=1 if Ij⊂fIi and aij=0 if IjfIi.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 fk (kth-iterate of f) is (17)sf=limk→∞lfkkand the topological entropy of a unimodal interval map f, denoted by Tentf, is given by (18)Tentf=log⁡sf=log⁡δmaxMf,where δmax(Mf) is the spectral radius of the transition matrix Mf ([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 (RLLLLC)∞. Putting the orbital points in order we obtain(19)x2<x3<x4<x5<x0<x1.The corresponding transition matrix is(20)Mf=0100000100000100000111111which has the characteristic polynomial (21)pt=det⁡Mf-tI=1+t+t2+t3+t4-t5.The growth number s(f) (the spectral radius of matrix M(f)) is 1.96595…. Therefore, the value of the topological entropy can be given by (22)Tentf=log⁡sf=0.675975….The study of the kneading sequences allows us to identify symbolic periodic orbits. In particular, with periodsn≤5, we have the following: 5-period - RLLLC∞, 4-period - RLLC∞, 5-period - RLLRC∞, 3-period - RLC∞, 5-period - RLRRC∞, 4-period - RLRC∞, 2-period - RC∞, and 1-period - C∞. 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.Table 1 Kneading sequences Characteristic polynomial Topological entropy RC 1 - t 0 RLRC - 1 + t + t 2 - t 3 0 RLRRC - 1 + t - t 2 - t 3 + t 4 0.414013 … RLC - 1 - t + t 2 0.481212 … RLLRC 1 - t - t 2 - t 3 + t 4 0.543535 … RLLC 1 + t + t 2 - t 3 0.609378 … RLLLC - 1 - t - t 2 - t 3 + t 4 0.656255 …The variation of the topological entropy in the parameter region depicted in Figure2(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 h. 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 h, 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,Sn—is associated with regimes of chaotic behavior.It is fundamental to emphasize that the existence of one-dimensional 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 maximumSnmaxof the densities of susceptible predicts the following pickSn+1maxof the same population.As pointed out before, the data from the chaotic time series (Figure2(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 Sn+1max=fSnmax allows us to tell a lot about the dynamics: given S0max, we canpredictS1max by S1max=fS0max, then use that information to predict S2max=fS1max, 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ε.


## 2.1. 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 h,Sn,In-bifurcation diagram and in Figure 1(b) the h,Sn-bifurcation diagram as a result of the variation of the time step h, giving a particular emphasis to the dynamical variable Sn. When the time step h is sufficiently small, the dynamical behaviors of the discrete-time model (5) are similar to the ones exhibited by the continuous-time SIR epidemic model. At increasing values of h, complex behaviors appear, namely, a period-doubling route to chaos (the variable In has a similar dominant qualitative behavior to Sn). The threshold value of the time step h 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 h goes beyond the central objectives of the present study.Figure 1
(a) 3D bifurcation diagram corresponding to the dynamical variablesSn and In, with the time step h as the control parameter. (b) Dynamics of the variable Sn, taking also h as the bifurcation parameter. (c) Variation of the maximal Lyapunov exponent of the system (5) with h. In all the previous situations β(t)=β0=0.445 and h∈2.92,3.4. The displayed region corresponds to the interval of values h, h∈M=mi,mf=3.25500,3.31097, and is associated with noticeable features of the dynamical behavior of Sn, as we will see in our study below.
(a) (b) (c)In the context of epidemic models, it is usual to define the so-calledbasic reproduction number, Rp. In the case of constant contact rate, β(t)=β0, it is given by (6)Rp=β0μ+γ,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) (7)S0∗,I0∗=αμ,0,corresponding to the population with no infected individuals. In this case, the epidemic cannot maintain itself in the population and the DFE point S0∗,I0∗ is locally stable.(ii) IfRp>1, then model (5) also has an endemic equilibrium (EE) (8)S1∗,I1∗=γ+μβ0,-μβ0+αγ+μ,corresponding to a self-sustaining group of infected individuals in the population. In this case, the EE point S1∗,I1∗ 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 h. 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–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:(9)yn+1=Fyn,where F:Rm→Rm is a differentiable vector function and y is a m vector. Let DFy denote the matrix ∂Fi/∂yi of partial derivatives of the components Fi at y. Then the corresponding matrix of partial derivatives for the nth iterate Fn of F is given by (10)∂Fni∂yj=DFFn-1y⋯DFFyDFy.Considering(11)An=DFFn-1y⋯DFFyDFy,we have(12)An=DFyn-1⋯DFy1DF0.Then(13)An+1=DFynAn.In this context,An can be decomposed into a product of an orthogonal matrix and an upper triangular matrix with positive diagonal elements. More precisely,(14)Qn+1Rn+1=DFynQnRn,where Qn, Qn+1, Rn, and Rn+1 are the upper triangular matrices with positive diagonal elements. Taking Rijn,i=1,…,m as the set of diagonal elements of Rn, the Lyapunov exponents are given by (15)λi=limn→∞ln⁡Rijnn,with i=1,…,m.We display in Figure1(c) the variation of the maximal Lyapunov exponent of system (5), λmax, with h 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 h in the interval 3.014,3.182, a phenomenon of homeochaos was found corresponding to 0<λmax≪1.The region in red displayed in Figure1(b) and in Figure 1(c), corresponding to values of the time step h in the interval M=mi,mf=3.25500,3.31097, is associated with salient, eye-catching, and noteworthy features of the dynamics of the variable Sn that will be studied in the following section.


## 2.2. Iterated Maps, Symbolic Dynamics, and Topological Entropy

In this paragraph, a considerable relevance is given to the dynamical variableSn studied for h within the interval M=mi,mf=3.25500,3.31097, presented previously. We show in Figure 2(a) a typical SI-phase diagram and in Figure 2(b) the corresponding variation of the iterates of Sn for h=3.31075. In order to gain additional significant qualitative insights into the principles and mechanisms underlying the dynamics of these iterates of Sn, we started our study by recording the successive relative (local) maxima of Sn.Figure 2
(a)SnIn-phase diagram (h=3.31075). (b) Series of the iterations of Sn. (c) Unimodal-type iterated map corresponding to the successive points Snmax,Sn+1max, where Snmax denotes the nth local maximum (h=3.31075). (d) Bifurcation diagram of Snmax, taking h∈M=3.25500,3.31097 as the control parameter. € Variation of the topological entropy corresponding to the dynamics of Snmax for h∈M=3.25500,3.31097.
(a) (b) (c) (d) (e)Each local maximum corresponds to a peak in the susceptible series of iterates. Considering the pairsSnmax,Sn+1max, where Snmax denotes the nth local maximum, an unimodal-type iterated map emerges. As shown in Figure 2(c) the data from the series of values of Sn, corresponding to successive points Snmax,Sn+1max, appear to fall on a logistic curve. Indeed, treating the graph as a representation of a function Sn+1max=fSnmax 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 I=a,b 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 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–31]).Let us consider a unimodal mapf and the interval I subdivided into the sets IL=[a,c[, IC={c}, and IR=]c,b]. The function f is monotonically increasing for x∈IL and monotonically decreasing for x∈IR (Figure 2(c)). We call a lap of f each of such maximal intervals where the map f is strictly increasing or strictly decreasing. The total number of distinct laps is called the lap number of f and it is usually denoted by l=lf. The left and the right subintervals are labeled by the letters L and R, respectively, and the set IC will be denoted by C. The symbolic sequence starting from f(c) 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 f, Oc, obtained by iterating the map, that is, Oc=xi:xi=fic,i∈N. From this numerical orbit, Oc, we get a symbolic sequence (symbolic orbit) S=S1S2…Sj…, where(16)Sj=Liffjx<c,Sj=Ciffjx=cSj=Riffjx>c,that characterizes the dynamics. When O(c) is a k-periodic orbit, we obtain a sequence of symbols that can be characterized by a block of length k, the kneading sequence S(k)=S1S2…Sk-1C. The orbit Oc, which is made of ordered points xi, determines a partition P(k-1) of the invariant range I=f2(c),f(c)=x2,x1 into a finite number of subintervals labeled by I1,I2,…,Ik-1. This partition is associated with a (k-1)×(k-1) transition matrix M=[aij], where the rows and columns are labeled by the subscript of subintervals and the matrix elements are defined as aij=1 if Ij⊂fIi and aij=0 if IjfIi.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 fk (kth-iterate of f) is (17)sf=limk→∞lfkkand the topological entropy of a unimodal interval map f, denoted by Tentf, is given by (18)Tentf=log⁡sf=log⁡δmaxMf,where δmax(Mf) is the spectral radius of the transition matrix Mf ([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 (RLLLLC)∞. Putting the orbital points in order we obtain(19)x2<x3<x4<x5<x0<x1.The corresponding transition matrix is(20)Mf=0100000100000100000111111which has the characteristic polynomial (21)pt=det⁡Mf-tI=1+t+t2+t3+t4-t5.The growth number s(f) (the spectral radius of matrix M(f)) is 1.96595…. Therefore, the value of the topological entropy can be given by (22)Tentf=log⁡sf=0.675975….The study of the kneading sequences allows us to identify symbolic periodic orbits. In particular, with periodsn≤5, we have the following: 5-period - RLLLC∞, 4-period - RLLC∞, 5-period - RLLRC∞, 3-period - RLC∞, 5-period - RLRRC∞, 4-period - RLRC∞, 2-period - RC∞, and 1-period - C∞. 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.Table 1 Kneading sequences Characteristic polynomial Topological entropy RC 1 - t 0 RLRC - 1 + t + t 2 - t 3 0 RLRRC - 1 + t - t 2 - t 3 + t 4 0.414013 … RLC - 1 - t + t 2 0.481212 … RLLRC 1 - t - t 2 - t 3 + t 4 0.543535 … RLLC 1 + t + t 2 - t 3 0.609378 … RLLLC - 1 - t - t 2 - t 3 + t 4 0.656255 …The variation of the topological entropy in the parameter region depicted in Figure2(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 h. 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 h, 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,Sn—is associated with regimes of chaotic behavior.It is fundamental to emphasize that the existence of one-dimensional 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 maximumSnmaxof the densities of susceptible predicts the following pickSn+1maxof the same population.As pointed out before, the data from the chaotic time series (Figure2(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 Sn+1max=fSnmax allows us to tell a lot about the dynamics: given S0max, we canpredictS1max by S1max=fS0max, then use that information to predict S2max=fS1max, 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ε.


## 3. 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 3D-plot that illustrates the variation of the maximal Lyapunov exponent within theε,h parameter space.To start with, we exhibit in Figure3 the graph of the step function β(t) 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 (h,Sn,In) bifurcation diagram and in Figure 4(b) the (h,Sn) bifurcation diagram, as a consequence of the variation of the time step h. 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 t. 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 h, giving rise to a particular set of maximal Lyapunov exponents, denoted by MLE. With the 3D graph of Figure 5, we intend to illustrate the variation of the maximal Lyapunov exponent in the (h,ε) parameter space. In this context, for each ε we have a particular set MLE, obtained with the variation of the time step h. 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 h∈3.14,3.17 and characterized by 0<λmax≪1.Figure 3
Graph ofβ(t) for ε=0.002, T=0.5, and t∈0,4.Figure 4
(a) 3D bifurcation diagram corresponding to the dynamical variablesSn and In, with the time step h as the control parameter. (b) Dynamics of the variable Sn, taking also h as the bifurcation parameter. (c) Variation of the maximal Lyapunov exponent of the system (5) with h. In all the previous situations ε=0.002 and h∈2.92,3.4.
(a) (b) (c)Figure 5
3D graph of the maximal Lyapunov exponent, with0<ε≤0.0075 and 2.92≤h≤3.4.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.


## 4. 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, andv0 is the proportion of those vaccinated successfully (with 0<v0<1). Indeed, constant vaccination takes action in the first equation of model (5), reducing the birth rate, α, of susceptible, so that the system becomes(23)Sn+1=Sn+h1-v0α-μSn-βnhSnInIn+1=In+hβnhSnIn-γ+μIn,where the term(1-v0)α=α-v0α incorporates the density of newborn infants vaccinated, v0α. 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+∈MLE+, with MLE+ the set of all values λmax+ that result from the variation of the time step h, for fixed values of ε and v0. 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, v0, for a lower degree of seasonality and also for a higher degree of seasonality ε. The respective derivatives of λ¯max+ in order to v0, dλ¯max+/dv0, are depicted in Figure 6(b). A prominent feature emerges from these two graphs of Figure 6:the value ofv0=0.02seems to act as a threshold, having a decisive effect on the dynamics. In this context, for values v0<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(a). For values v0≥0.02 those types of points tend to disappear (please see the circles of Figures 7(b) and 7(c)). As a consequence, with the chaotic behavior suppressed, the system converges to controlled population densities. This threshold property of v0 is illustrated in Figure 7 taking, as an example, ε=0.002. All our numerical simulations for other values of ε also confirm v0=0.02 as a threshold of the dynamics (results not shown).Figure 6
(a) Variation of the mean value of the positive Lyapunov exponents,λ¯max+, with the proportion of susceptibles who received constant vaccination, v0, for ε=0.001 (green) and ε=0.0055 (blue). (b) Variation of the derivative function dλ¯max+/dv0, for ε=0.001 (green) and ε=0.0055 (blue).
(a) (b)Figure 7
Behavior of the maximal Lyapunov exponentλmax for different values of the constant vaccination v0: (a) v0=0.015; (b) v0=0.02 (thethresholdvalue of v0); and (c) v0=0.02.
(a) (b) (c)


## 5. 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 steph 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 unimodal-type 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,Sn. 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.



---

*Source: 7191487-2018-12-02.xml*