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.
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 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 equationsThe 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) takeswhere gives the mean contact rate (the base of transmission rate). All the parameters of system (1) that have been described, i.e., , , , and , are nonnegative and smaller than 1. The parameter , with , 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 asIn this context, we have in the high season and 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 modelwhere 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 modelwhich 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, . 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 , taking , (Section 2) and epidemics with seasonal contact rate , taking , with given by (3) (Section 3).
In our analysis, if not otherwise specified, we will set , , , . A variety of interesting dynamical behaviors will occur for values of the strength of the seasonal forcing, , such that and for the time step as a control parameter, with , 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, , 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 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 and , 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 . 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 -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 discrete-time 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.
(a)
(b)
(c)
In the context of epidemic models, it is usual to define the so-called basic reproduction number, . In the case of constant contact rate, , 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 is locally stable.
(ii) If , 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 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–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 is a differentiable vector function and is a vector. Let denote the matrix of partial derivatives of the components at . Then the corresponding matrix of partial derivatives for the iterate of 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 , , , and are the upper triangular matrices with positive diagonal elements. Taking as the set of diagonal elements of , the Lyapunov exponents are given by with
We display in Figure 1(c) the variation of the maximal Lyapunov exponent of system (5), , with along the interval . In agreement with the bifurcation diagrams, we have within the periodic windows and 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 , a phenomenon of homeochaos was found corresponding to .
The region in red displayed in Figure 1(b) and in Figure 1(c), corresponding to values of the time step in the interval , is associated with salient, eye-catching, and noteworthy features of the dynamics of the variable 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 variable studied for within the interval , 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 . 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 .
(a)
(b)
(c)
(d)
(e)
Each local maximum corresponds to a peak in the susceptible series of iterates. Considering the pairs , where denotes the local maximum, an unimodal-type iterated map emerges. As shown in Figure 2(c) the data from the series of values of , corresponding to successive points , appear to fall on a logistic curve. Indeed, treating the graph as a representation of a function 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 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 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, From this numerical orbit, , we get a symbolic sequence (symbolic orbit) , wherethat 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 . The orbit , which is made of ordered points , determines a partition of the invariant range into a finite number of subintervals labeled by . This partition is associated with a transition matrix , where the rows and columns are labeled by the subscript of subintervals and the matrix elements are defined as if and if .
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 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 obtainThe corresponding transition matrix iswhich has the characteristic polynomial The growth number (the spectral radius of matrix ) is . 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 , we have the following: 5-period - , 4-period - , 5-period - , 3-period - , 5-period - , 4-period - , 2-period - , 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 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 maximum of the densities of susceptible predicts the following pick 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 allows us to tell a lot about the dynamics: given , we can predict by , then use that information to predict , 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 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 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 , 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, , 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 (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 and characterized by .
(a)
(b)
(c)
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, and is the proportion of those vaccinated successfully (with ). 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 incorporates the density of newborn infants vaccinated, . 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 , with the set of all values that result from the variation of the time step , for fixed values of and . In Figure 6(a) we represent the behavior of the mean value of the positive Lyapunov exponents, , as a function of the proportion of susceptibles who received constant vaccination, , for a lower degree of seasonality and also for a higher degree of seasonality . The respective derivatives of in order to , , are depicted in Figure 6(b). A prominent feature emerges from these two graphs of Figure 6: the value of seems to act as a threshold, having a decisive effect on the dynamics. In this context, for values the dynamics is still exhibiting chaotic behavior, with for some time steps corresponding to the points inside the dashed circle of Figure 7(a). For values 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 is illustrated in Figure 7 taking, as an example, . All our numerical simulations for other values of also confirm as a threshold of the dynamics (results not shown).
(a)
(b)
(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 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 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, . 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.
Data Availability
The research data used to support the findings of this study are included within the article. Our paper is totally self-contained, with clearly cited supporting references from the literature. In the context of our study, we provide an approach that we believe is still lacking in the literature. This analysis is supported by numerical results (obtained using Wolfram Mathematica 11.1.1) that match the theoretical predictions. The manuscript is technically correct and it is sufficiently accessible to the nonspecialist reader. There is a logical flow of argument and validity of conclusions drawn. This study provides another illustration of how an integrated approach, involving numerical evidence and theoretical reasoning, can directly enhance our understanding of biologically motivated models. It will probably stimulate the development of new research on epidemiology.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was partially supported by the Portuguese Foundation for Science and Technology (FCT). Particularly, within the projects UID/MAT/04106/2013 (CIDMA) (Cristina Januário) and UID/MAT/04459/2013 (Jorge Duarte and Nuno Martins).