Controlling infectious diseases: the decisive phase effect on a seasonal vaccination strategy
Jorge Duarte, Cristina Januário, Nuno Martins, Jesús Seoane, Miguel AF Sanjuán
aa r X i v : . [ m a t h . D S ] F e b Controlling infectious diseases: the decisivephase effect on a seasonal vaccination strategy
J. Duarte , ∗ , C. Janu´ario , N. Martins , J. M. Seoane , M. A. F. Sanjuan ISEL-Engineering Superior Institute of Lisbon, Department of Mathematics,Rua Conselheiro Em´ıdio Navarro 1, 1950-007 Lisboa and Center for MathematicalAnalysis, Geometry and Dynamical Systems, Instituto Superior T´ecnico,Universidade de Lisboa, Portugal Center for Research and Development in Mathematics and Applications(CIDMA) Department of Mathematics, University of Aveiro, 3810-193 Aveiro,Portugal Department of Mathematics and Center for Mathematical Analysis, Geometryand Dynamical Systems, Instituto Superior T´ecnico, Universidade de Lisboa, Av.Rovisco Pais 1, 1049-001 Lisboa, Portugal Nonlinear Dynamics and Chaos Group, Departamento de F´ısica, UniversidadRey Juan Carlos, Tulip´an s/n 28933 M´ostoles, Madrid, Spain ∗ E-mail: [email protected] (corresponding author)
Abstract
The study of epidemiological systems has generated deep interest in exploring thedynamical complexity of common infectious diseases driven by seasonally varyingcontact rates. Mathematical modeling and field observations have shown that, underseasonal variation, the incidence rates of some endemic infectious diseases fluctuatedramatically and the dynamics is often characterized by chaotic oscillations in theabsence of specific vaccination programs. In fact, the existence of chaotic behaviorhas been precisely stated in the literature as a noticeable feature in the dynam-ics of the classical Susceptible-Infected-Recovered (SIR) seasonally forced epidemicmodel. However, in the context of epidemiology, chaos is often regarded as an unde-sirable phenomenon associated with the unpredictability of infectious diseases. Asa consequence, the problem of converting chaotic motions into regular motions be-comes particularly relevant. In this article, we consider the phase control techniqueapplied to the seasonally forced SIR epidemic model to suppress chaos. Interestingly,this method of controlling chaos takes on a clear meaning as a weak perturbationon a seasonal vaccination strategy. Numerical simulations show that the phase dif-ference between the two periodic forces - contact rate and vaccination - plays a veryimportant role in controlling chaos.
Key words:
Infectious diseases, Seasonally forced SIR model, Vaccinationstrategy, Phase control
Preprint submitted to Elsevier Science 17 February 2021
Introduction
The outbreaks and spread of infectious diseases are a threat to public healthand a source of serious problems to the economic and social development of oursociety. The research in epidemic dynamics is, without doubt, critical to anyattempt to prevent or minimize the transmission of diseases. In particular, wehave all witnessed a remarkable intensification of this work in the present daysgiven the quantity of infected people worldwide due to the dissemination of anew virus, referred to as COVID-19 (Coronavirus disease 2019), an infectiousdisease emerged from China in November 2019 [2].The use of mathematical models has contributed greatly to our better under-standing of the underlying mechanisms that influence the spread of diseasesand has suggested control strategies, which may have significant managementimplications [3]. More precisely, these models are truly significant in differentfields, such as policy making, risk assessment, emergency planning and defi-nition of health-economic control-programs. The mechanism of transmissionof infections is now known for most diseases. Generally, diseases transmittedby viral agents (such as measles, influenza, rubella and chickenpox) conferimmunity against reinfection, while diseases transmitted by bacteria (such asmeningitis and tuberculosis) confer no immunity against reinfection. Otherdiseases, such as malaria, are transmitted not directly from human to humanbut by vectors (usually insects). The goal of mathematical epidemiology hasbeen the development of mathematical models for the spread of disease as wellas tools for their analysis.The first mathematical model describing an infectious disease was proposedby Bernoulli, in 1760, to study the spread of smallpox, which was prevalentat the time [4]. An early recognition of the importance of mathematical mod-eling came in 1911 when Ronald Ross won the Nobel Prize in Physiologyor Medicine, for his demonstration of the dynamics of the transmission ofmalaria between mosquitoes and humans (using models formulated with dif-ferential equations) [5]. Recent studies in epidemiology are particularly focusedon modeling the previously mentioned COVID-19 and high impact recurrentepidemics, best exemplified by childhood infectious diseases such as measles,chickenpox, mumps, whooping cough and rubella, but also including hepatitis,different types of influenza (see [6] and references therein).One of the questions that first attracts the attention of scientists interestedin the study of the spread of communicable diseases is how severe will an epi-demic be. This question may be interpreted in a variety of ways and studiedwith the aid of models. To formulate dynamic models predicting the behaviorof outbreaks and the transmission of infectious diseases, compartmental mod-els are usually considered, with the population in a given region subdivided2nto several distinct groups or compartments, depending upon their experiencewith respect to disease, whose sizes change with time. These compartments,which are mutually exclusive categories based on infection status, were ini-tially proposed by William O. Kermack and Anderson G. McKendrick in 1927with their mathematical epidemic nonlinear system of differential equations,called the Susceptible-Infected-Recovered (SIR) model [7]. In the context ofthis prominent model, of great historical importance in research of epidemics,the population being studied is labeled into three classes: the susceptible com-partment refers to individuals who have never come into contact with the dis-ease at time t , but they could catch it, i.e., they are vulnerable to exposurewith infectious people. Infected individuals, assumed infectious and capableof spreading the disease to those in the susceptible category, and remain inthe infectious compartment until their recovery. The recovered class refers toindividuals who have been infected and then recovered from the disease. Therecovered individuals are immune for life and are not able to transmit the in-fection to others. They are essentially removed from the population and playno further role in the dynamics. These models can include, among other as-pects, time-dependent parameters to represent the effects of seasonality andhuman demographics by adopting birth and death rates. In practice, given amodel, its parameters must be determined to represent a particular epidemiccontext. In formulating models as systems of differential equations, we areassuming that the epidemic process is deterministic, that is, the behavior ofa population is determined completely by its history and by the rules whichdescribe the model. Particularly, in the present study, we will consider a SIRmodel representing a class of seasonally forced epidemic models with vitaldynamics (birth and death rates) and constant population.The article is organized as follows. Immediately after the present Section I,where an introduction to our study is provided, Section II gives essential pre-liminaries about the model, as well as the control procedure used. Section IIIis devoted to a diagnosis of the existence of chaos in the system. We exhibitdifferent dynamical interactions between the densities of the susceptible andthe infected individuals varying the degree (or amplitude) of seasonality, ε .Revealing insights about the chaotic dynamics of the SIR model, incorporat-ing seasonal fluctuation, are gained through the computation of bifurcationdiagrams and the largest Lyapunov exponent. In Section IV, the phase controltechnique is applied in order to suppress chaos. A clear biological meaning isassigned to this method of controlling chaos as a weak/small perturbation ona seasonal vaccination strategy. Finally, our last considerations are devoted tosignificative conclusions. 3 Description of the model
In the literature, and as far as the description of the model is concerned, tworelated sets of time-dependent variables have been considered in the epidemicmodeling process. The first set of dependent variables counts the number ofpeople in each of the groups, each as a function of time. Considering that thepopulation size ( N ) is defined as N = e S + e I + e R , the dynamical behaviorof a homogeneously mixing population is described using the compartmentalforced e S e I e R model [8], with vital dynamics (birth and death rates), given by d e Sdt = σ − µ e S − β ( t ) e S e INd e Idt = β ( t ) e S e IN − ( γ + µ ) e I d e Rdt = γ e I − µ e R . (1)The second set of dependent variables represent a fraction of the total popula-tion in each of the three categories. So, being N the previous total population,we have S = e S ( t ) N , the susceptible fraction of the population, (2) I = e I ( t ) N , the infected fraction of the population, R = e R ( t ) N , the recovered fraction of the population,verifying S + I + I = 1. The two sets of dependent variables are proportionalto each other, giving the same information about the progress of the epidemic.The independent variable is time t , measured in a specific unit according tothe epidemic context (such as days, years, etc.). Although it may seem morenatural to work with the populations counts, some of our calculations willbe simpler if we use the fractions instead. As a consequence, we are goingto consider the SIR model with vital dynamics and constant population ( ∀ t, dSdt + dIdt + dRdt = 0), such that S + I + R = 1 , dSdt = σ − µS − β ( t ) SINdIdt = β ( t ) SIN − ( γ + µ ) I dRdt = γI − µR , (3)with time t scaled in unit of years. The parameters of the model (3) andrespective meaning are presented in the following table (for more details, pleasesee [3], [9], [10] and references therein).4able 1. Description of the parameters of model (3)Parameter (year) − Description σ Birth rate µ Natural death rate γ Recovery rate β ( t ) Contact or transmission rateWe assume that the newborns are susceptible and that the birth and deathrates are balanced, σ = µ . As we will exemplify with particular values ofthe parameters, the average time to recover from infection can be derivedfrom γ and given by 1 /γ . The contact or transmission rate of the infection(the coefficient of infectivity) β ( t ) represents the number of contacts withother individuals per infective per unit of time. More intensive research onseasonally-forced epidemic models did not begin until the 1970s ([11], [12]).With the purpose of expressing seasonality, a commonly used scheme for thecontact rate β ( t ) takes the form β ( t ) = β (1 + εϕ ( t )) , (4)where β gives the mean contact rate; parameter ε , with 0 ≤ ε ≤
1, representsthe strength 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.Modelers often incorporate seasonality by making the contact rate β ( t ) a si-nusoidal function of time β cos ( t ) = β (1 + ε cos (2 πt )). With this procedure,a seasonally-varying transmission rate yields oscillations at periods that areinteger multiples of the period of forcing.As recently proposed in [13], and inspired by the developed arguments, weare going to consider in our epidemiological study yearly periodic Kot-typefunctions. Particularly, the previously mentioned function β ( t ) is given by β ( t ) = β ε + cos (2 πt )1 + cos (2 πt ) !! . (5)In order to point out the main features of the adopted type of periodic func-tions, we compare in Fig. 1 a sinusoidal-type periodic function with the pe-riodic Kot-type function β ( t ). Differently from the periodic sinusoidal-typefunction, the Kot-type function gives asymmetrical weights to the seasonalregimes, stressing the relative maxima of the contact rate. This is an eye-catching and noteworthy feature that has been pointed out as more realistic.An exposure of a susceptible to an infectious is an encounter in which the5 β ( t ) , β c o s ( t ) Fig. 1. Comparison between the periodic Kot-type function β ( t ) (Black) and theusual periodic function β cos ( t ) (Blue) ( ε = 0 . infection is transmitted. In this context, the contact rate, β ( t ), is the averagedensity of susceptible in a given population contacted per infectious individualsper unit of time. Therefore, β ( t ) S ( t ) denotes the rate of the total density ofsusceptible infected by one infectious and β ( t ) S ( t ) I ( t ) represents the rate ofinfection of the susceptible by all infectious.It is well known that seasonal forces, including climatic factors and humanphenomena (such as school schedules), are a primary factor responsible forthe transmission and dynamical behavior of most recurrent infectious diseases[14,15,16,17,18,19,20]. As a consequence, it becomes natural to model thesediseases as periodically forced nonlinear systems ([18,21,22,23,24]). Seasonalforces in these systems shape the spread of the infectious disease and manystudies have stated that intense seasonality can lead to strong erratic patternswith the presence of chaotic oscillations ([1] and [22]). The chaotic behavior ofthe simple SIR model with seasonality has deep biological consequences. Thesensitive dependence on the initial conditions, in the final epidemic outcome,is one of the main concerns. It is important to notice that this rich dynamicalscenario is a direct consequence of seasonality, which makes any vaccinationstrategy difficult to design. In fact, the inclusion of seasonality increases thecomplexity of the models at several levels, making a complete comprehensionof its influence a non-trivial and challenging task. In this context, chaos isregarded as an unwanted feature and, as a consequence, control schemes (con-sidered as methods of controlling chaos, suppressing chaos, or taming chaos)play a prominent/decisive role. In any case, the explicit aim of all the pro-cedures is to obtain stable periodic orbits from chaotic ones by applying asmall/weak, and carefully chosen, perturbation to the system.The controlling methods have been traditionally classified into two generalcategories: feedback methods and nonfeedback methods, based on their inter-action with the system’s dynamics. Feedback methods attempt to suppresschaos by stabilizing orbits already existing in the chaotic attractor. Nonfeed-back methods have been essentially used to stabilize periodically driven chaoticdynamical systems. Typically, they make stable solutions to appear by apply-ing small driving forces directly to some of the parameters of the system oras an additional forcing. It has been shown in the literature (for instance, in6apers [25,26,27,28,29], and references therein) that a phase difference φ , be-tween the main driving and the perturbation, influences profoundly the globaldynamics of the system.As we will see in the following lines, given the nature of the forced SIR epi-demic model, it is specifically a nonfeedback procedure, that makes use of theproperty where φ acts as a control parameter, that stands out to be partic-ularly useful in controlling chaotic behavior. Acting as a positive/beneficialbiological consequence, this nonfeedback scheme called phase control of chaosmakes a vaccination strategy clearly easier to design.Following the previous considerations, the particular system we are going toconsider in this article is dSdt = σ − µS − β ( t ) SI − v ( t ) S dIdt = β ( t ) SI − ( γ + µ ) I dRdt = γI − µR + v ( t ) S , (6)where the first periodic force, the contact rate β ( t ), drives the system to thechaotic state, while the second one, the vaccination rate v ( t ), is a weak periodicperturbation sensitively modifying the system’s dynamics. We assume that thevaccination rate v ( t ) is defined by the Kot-type function v ( t ) = v + α + cos (2 πrt + φ )1 + cos (2 πrt + φ ) ! , (7)where v gives the mean vaccination rate, φ is the phase difference between theapplied perturbation v ( t ) and the driving force β ( t ), a key parameter for ourcontrol scheme. Parameter r is the ratio of the frequency of those forces and α ,with α <<
1, measures the degree (or amplitude) of the seasonality of v ( t ).Once vaccinated, the individuals are no more susceptible (fact representedby the term − v ( t ) S in the first equation of the model) and the vaccinatedindividuals are added to the recovered class (fact represented by the term+ v ( t ) S in the third equation of the model).The vaccination switches individuals directly from the susceptible state ( S )to the immune state ( R ). We suppose that infection transmission rate β ( t )is subjected to an environmental periodic forcing and management imposesa countercycle control strategy through v ( t ). The rationale underpinning avaccination policy is to ensure that the proportion of susceptible individualsin the population would stay below a certain threshold. We display in Fig. 2 ajoint representation of the yearly periodic functions β ( t ) and v ( t ). Notice that,for φ >
0, the qualitative behavior of v ( t ) anticipates the qualitative behaviorof β ( t ) . In Section III, we will see how this behavior influences drastically the7ynamics, acting as a controller to suppress chaos. Figure 2 also suggests thatthis decisive effect of the phase difference will be achieved with just a smallperturbation, an observable fact reflected by the scale differences between β ( t )and v ( t ) ( v ( t ) << β ( t )). - - - (cid:1) ( t ) v ( t ) Fig. 2. Comparison between the contact rate Kot-type function β ( t ) (Black) and thevaccination Kot-type function v ( t ) (Red) ( ε = 0 . α = 0 . r = 2 and φ = π ). The model is clearly understood when β ( t ) = β and v ( t ) = 0. For this caseof constant contact rate, and in the absence of a vaccination strategy, thedynamical system (6) has two equilibrium points:(i) The Disease-Free Equilibrium (DFE)( S ∗ , I ∗ , R ∗ ) = σµ , , ! , (8)corresponding to a population with no infected individuals;(ii) The Endemic Equilibrium (EE)( S ∗ , I ∗ , R ∗ ) = γ + µβ , µβ ( R − , γβ ( R − ! , (9)corresponding to the case in which there is a significant group of infectiousindividuals, where R = βσµ ( µ + γ ) i.e., R = β ( µ + γ ) , for σ = µ ! (10)is a derived basic reproduction number (or basic reproductive number) withthreshold properties. The basic reproduction number of an infection, R , is ameasure of the potential for the disease to spread in a population. Heuristi-cally, R may be read as the expected number of secondary infections directlygenerated by a single infectious individual in a wholly susceptible population.The number R is not a biological constant for a pathogen as it is also af-fected by other factors such as the behavior of the infected population andenvironmental conditions. It can also be modified by physical distancing and8ther public policy or social interventions. R values are usually estimatedfrom mathematical models, and the estimated values are dependent on themodel used and values of other parameters. Therefore, values given in theliterature only make sense in the given context and it is recommended not tocompare directly values based on different models. The most important usesof R are determining if an emerging infectious disease can spread in a popula-tion and determining what proportion of the population should be immunizedthrough vaccination in order to eradicate a disease. Generally speaking, andindependently from biologically meaningful initial values, this means that:(a) If R <
1, the epidemic cannot maintain itself, since each infected individ-ual on average infects less than one member of the population. The EE point( S ∗ , I ∗ , R ∗ ) is then unstable, while DFE ( S ∗ , I ∗ , R ∗ ) is locally stable and thedisease goes extinct;(b) If R >
1, each infected individual infects more than one other member ofthe population and a self-sustaining group of infectious individuals will prop-agate. In this case, the EE point ( S ∗ , I ∗ , R ∗ ) is locally stable, while the DFEpoint ( S ∗ , I ∗ , R ∗ ) is unstable. The disease will remain permanently endemicin the population.Mathematical modeling in epidemiology provides understanding of the under-lying mechanisms that influence the spread of the disease and, in the process,it suggests control strategies. One of the significant results in mathematicalepidemiology is that the majority of mathematical epidemic models, includingthose that have a high degree of heterogeneity, usually exhibit threshold be-havior. In epidemiological terms, this can be stated as follows: If the averagenumber of secondary infections caused by an average infectious, the previouslycalled basic reproduction number, is less than one a disease will die out, whileif it exceeds one there will be an epidemic. In the context of differential equa-tion models (or, more generally, evolution equation models), R arises as adimensionless number of transmission. Throughout this work, we will use ourparameters following previous studies in the literature (see [3], [9], [10] andreferences therein). As a prior notice, it is important to emphasize that, in theframework of our study, the main qualitative features of the dynamics are notsensitive to the precise values of certain parameters. With the time measuredin units of years, the values of the parameters, corresponding to an infectiousdisease, are listed in the following table.9able 2. List of the parameter values. σ = µ = 0 .
01 (year) − = ⇒ µ = 100 years(mean lifetime of the host) γ = 50 (year) − = ⇒ γ = 0 .
02 year ≈ β = 1505 (year) − = ⇒ ≈ − (mean number of contacts per day) R ≈ β ( µ + γ ) ≈ ε , of the mainforce β ( t ), as well as parameters v , α and φ of the vaccination rate (theadded small perturbation that will act as a control component), will be takenas control parameters.At this moment, it is also critical to stress that, throughout our study, a closeattention will be focused, not on the phase control method per se , but on itsbiological significance and on its consequences in an epidemiological context. The complexity of the dynamics increases considerably with the introductionof the seasonal contact rate. In this section, we will examine the long-termbehavior of the three-dimensional chaotic attractors arising in the forced SIRsystem (6), with a seasonal transmission rate component and under a constantvaccination strategy. The time series, displayed in Fig. 3 (Lower panel), governthe population dynamics of the susceptible and infected. They were obtainedusing two different, but close, initial conditions.In this paragraph, the Lyapunov exponents of the SIR model (6) with the sea-sonal component β ( t ) and v ( t ) = v , receive our attention as a framework todiagnose chaos in the system. A discussion about the Lyapunov exponents as aquantitative measure of the rate of separation of infinitesimally close trajecto-ries, as well as a computation method, can be found in [30]. The characteristicLyapunov exponents measure the typical rate of the exponential divergence ofnearby trajectories in phase space, i.e., they give us information on the rateof growth of a very small error on the initial state of the system.10 lready chaotic ϵ S Already chaotic ϵ I λ ≈ ( Chaos ) ϵλ t I Fig. 3. Global noticeable features of the dynamics of the SIR system (6) considering aconstant vaccination strategy, v ( t ) = v = 0 .
071 and α = 0. Increasing values of thedegree of seasonality of the contact rate, ε , bring the dynamics to a chaotic regime.Upper panel: bifurcation diagrams of the dynamical variables S and I , taking ε as control parameter (0 . ≤ ε ≤ . λ , with ε (0 . ≤ ε ≤ . S and I , obtained using two different (butclose) initial conditions displayed in red and green ( ε = 0 . A positive Lyapunov exponent is commonly taken as an indicator of chaoticbehavior. In Fig. 3, taking ε as a control parameter, we present the variationof the largest Lyapunov exponent λ (Middle panel) and two bifurcation dia-grams (Upper panel) regarding the dynamical variables S and I . In agreementwith the represented bifurcation diagrams, the largest Lyapunov exponent, λ ,is positive in the chaotic regime and λ ≈ ε results in the complexity of the dynamics at highervalues. This way, the actual existence of chaos is numerically recognized, withthe seasonal transmission function β ( t ) as the force driving the system to thechaotic state. In particular, for the choice of parameters β = 1505, ε = 0 . β ( t ) and v ( t ) = v = 0 . λ ≈ .
1. From now on, we keep the parametervalues ( β = 1505, ε = 0 .
138 and v = 0 . ϕ / πα r = ϕ / πα r = λ > < λ ⩽ (cid:6) λ ⩽ (cid:7) λ ⩽ N(cid:8)(cid:9) - c(cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15) r(cid:16)(cid:17)(cid:18)(cid:19)(cid:20) Fig. 4. Density plots for the largest Lyapunov exponent corresponding to r = 2, inthe ( φ, α )-parameter region. Left: 0 ≤ φ ≤ π and 0 . ≤ α ≤ .
01; Right: zoomin of the previous density plot considering π ≤ φ ≤ π and 0 . ≤ α ≤ . λ > As we have just stated, intense seasonality of the periodically varying contactrate induces chaotic dynamics in the epidemic system under a conventionalconstant vaccination strategy. However, in the context of epidemiology, chaosis often regarded as an undesirable phenomenon associated with the erraticpermanence of infectious diseases. As a consequence, the aim of the presentsection is to provide a comprehensive study of the forced epidemic system interms of the implementation of the phase control, as a control technique tosuppress the chaotic behavior.Inspired by the importance of vaccination for the elimination of infectiousdiseases (please see [31]), this procedure is introduced as a biologically mean-ingful weak perturbation on a seasonal vaccination strategy. More specifically,with the parameters β = 1505, ε = 0 .
138 and v = 0 .
071 already tailored insuch a way that the asymptotic state of system (6), under a constant vaccina-tion strategy ( v ( t ) = v ), is chaotic, our aim here is to analyze the effects ofthe phase φ in the chaotic regime, when an additive perturbation to v , givenby (7), is included. For this purpose, the amplitude of the perturbation α isassumed to be very small, i.e., α << φ , combined with the variation of α , we have displayed in Fig. 4 density plots. The color scale corresponds tothe computation of the largest Lyapunov exponent over the parameter regioncharacterized by 0 ≤ φ ≤ π and 0 . ≤ α ≤ .
01, fixing r = 2. In Fig. 4, wedistinguish two dynamical regimes. The chaotic regime, corresponding to λ >
0, represented by the colored region (Blue: 0 < λ ≤ . . < ϕ /π S 1 1.2 1.4 1.6 1.8 20.000240.000280.00032 ϕ /π I ϕ /π λ (cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26) (cid:27) I Fig. 5. Upper panel: bifurcation diagrams for S and I when α = 0 . φ ascontrol parameter ( π ≤ φ ≤ π ); Middle panel: variation of the largest Lyapunovexponent, with α = 0 .
009 and π ≤ φ ≤ π ; Lower panel: time series for the dy-namical variables S and I , for α = 0 .
009 and φ = π (a couple of parameter valuesmarks the beginning of a transition from chaotic to regular behavior). λ ≤ .
01; Re ˙d: 0 . < λ ≤ . λ > . r = 2, we are able to obtain aclear transition from chaotic to regular behavior increasing φ . In a biologicalcontext, this choice of r guarantees a frequency of the vaccination withinone year. With these density plots, we can better appreciate the structureof the chaotic region (colored) and periodic region (white). In our numericalexploration, the choice of such small values for the perturbation parameter α allows us to demonstrate the effectiveness of the phase control. The densityplot (right-hand side of Fig. 4), is precisely a zoom in of the previous densityplot (left-hand side of Fig. 4), which clearly exhibits the suppression of chaos,increasing the values of φ , when π ≤ φ ≤ π . In particular, the strikingcharacteristics of the phase control - the key role of the phase φ in selectingthe final state of the system successfully, combined with tiny values of α tosuppress chaos - are present in this density plot (right-hand side of Fig. 4). Inorder to continue illustrating this primary role of the phase φ we show, in Fig. 5(upper panel), bifurcation diagrams for S and I , fixing α = 0 .
009 and taking φ as control parameter, with π ≤ φ ≤ π and the corresponding variationof the largest Lyapunov exponent (Middle panel). In fact, we have found awide range of phase values producing regular motion, in which the systemleaves the chaotic region to periodic motion via inverse period doubling. A13emarkable feature is that for values of α , such that α % . φ (density plot of Fig. 4, right-hand side). In particular, for α =0 . φ ≈ π . This observation indicatesthat the phase φ is a sensitive parameter for the system bifurcation, i.e.,the distribution of regular and chaotic regions strongly depends on the phasedifference (please, visit again Fig. 2, and see the representation of v ( t ) preciselyfor φ ≈ π . It is extremely interesting that the weak additive force makes thevaccination strategy seasonal and its phase term (with φ >
0) ensures thatthe vaccination v ( t ) anticipates the seasonality of the transmission rate β ( t ).Therefore, the phase control method for suppressing chaos takes on a clearmeaning, where the phase difference between the two central driven forces:(i) a transmission rate and (ii) a vaccination strategy, stand out to be a keyperturbation parameter with immediate and beneficial biological consequences- the control of a given infectious disease. For decades, mathematical models have been proposed to evaluate the spreadand control of infectious diseases. Recently, for the last entire year, the worldhas been experiencing the intensification of the research in epidemic dynamicsgenerated by the dissemination of COVID-19 ([32],[33],[34],[35] and [36]). Theyear 2020 has seen significant advances taking place, to build the infrastruc-tures to keep up with this coronavirus. SARS-CoV-2, the virus which causesCOVID-19, is constantly evolving and mutating, as do all similar virus. Thepressure on the virus to evolve is increased by the fact that so many millionsof people have now been infected. Such identified changes since it emergedin 2019, in a RNA virus that exists as a cloud of genetic variants known asquasispecies, are completely to be expected to occur and have been useful inunderstanding the worldwide spread as well as the transmission patterns. Themajority of the mutations will not be significant or cause for concern, butsome may give the virus an evolutionary advantage which may lead to highertransmission. More specifically, a new genetic variant of the virus has emergedand is spreading in many parts of the UK and across the world. This efficienttransmission among people is usually associated to a modeling process with ahigher basic reproduction number, R .It is precisely in this context of controlling infectious diseases with R >> perse , but to its biological significance and to its consequences within the frame-14ork of epidemiology. We have provided detailed/revealing insights about therole played by the phase difference of the two periodically driven forces - theseasonal transmission function β ( t ) and a vaccination component v ( t ). Hav-ing started with the analysis of the existence of chaos, the transmission ratefunction β ( t ) has been identified as the main force driving the system to thechaotic regime under a classical constant vaccination strategy. With the com-putation of the largest Lyapunov exponent, values of the parameters of β ( t )have been tailored within the chaotic region.Motivated by the fact that, in an epidemiological context, chaotic behavioris often associated with the erratic permanence of infectious diseases and avaccination strategy is associated with their efficient elimination. Thus, an ideaemerged - to introduce the control procedure as a biologically meaningful weakperturbation on a seasonal vaccination function v ( t ). Given the importanceof controlling the chaotic behavior, i.e., the necessity of having predictabledensities for the epidemic populations, we have applied the mentioned periodiccontrol signal v ( t ) , including the phase difference with respect to the periodicforcing of the initial system, which has acted as an effective control strategy.More precisely, the chaotic epidemic outbreaks, that appeared as a result ofthe seasonal variations in the contact rate, have been suppressed by the usedvaccination control scheme. Indeed, the crucial role of the phase term, in theseasonal component of v ( t ), was evidenced by using numerical simulations thatallowed us to clearly identify dynamical transitions from a chaotic regime toa regular behavior, which is biologically associated with the control of a giveninfectious disease.This study provides another illustration of how an integrated approach, in-volving numerical evidences and theoretical reasoning, within the theory ofdynamical systems, can contribute to our understanding of important bio-logical models and provide a trustworthy explanation of complex phenomenawitnessed in biological systems.Above all, the recent appearance of COVID-19 is a reminder that there isstill so much to learn about a pandemic dynamics. The pace of the researcheffort in the past year has been extraordinary. However, there is no room forcomplacency. We have to be humble and mostly be prepared to adapt andrespond to new and continued changes. Acknowledgments
This research has been financially supported by the Portuguese Foundationfor Science and Technology (FCT) under Projects No. UIDB/04106/2020 andNo. UIDP/04106/2020 (CIDMA) (CJ) and No. UID/MAT/04459/2013 (JDand NM); and by the Spanish State Research Agency (AEI) and the EuropeanRegional Development Fund (ERDF, EU) under Projects No. FIS2016-76883-15 and No. PID2019-105554GB-I00 (JMS and MAFS).
References [1] Barrientos PG, Rodr´ıguez JA, Ruiz-Herrera A. Chaotic dynamics in theseasonally forced SIR epidemic model. J Math Biol 2017; 75:6-7. DOI10.1007/s00285-017-1130-9.[2] Libotte GB, Lobato FS, Platt GM, Neto AJS. Determination of anOptimal Control Strategy for Vaccine Administration in COVID-19 PandemicTreatment. Computer Methods and Programs in Biomedicine 2020; 196:105664.DOI: 10.1016/j.cmpb.2020.105664[3] Zhang Y, Zhang Q. Chaos analysis and control for a class of SIR epidemicmodel with seasonal fluctuation. Int J Biomath 2013; 6,1:1250063-1-11.https://doi.org/10.1142/S1793524512500635[4] Bernoulli D. Essai d’une nouvelle analyse de la mortalit´e caus´ee par la petiteverole et des avantages de l’inoculation pour la prevenir. Mem. Math. Phys.Acad. Roy. Sci. Paris 1760;1[5] Ross R, Howard LO, Gorgas WC. The prevention of malaria. 2 nd ed. London:John Murray; 1911.[6] Olinky R, Huppert A, Stone L. Seasonal dynamics and thresholds governingrecurrent epidemics. J Math Biol 2008; 56:827-839.[7] Kermack WO, Mckendrick AG. A contribution to the mathematical theory ofepidemics. Proc Roy Soc London A 1927; 115:700-721.[8] Dietz K. The incidence of infectious diseases under the influence of seasonalfluctuations. In: Mathematical models in medicine, Berlin: Springer; 1976, p.1-15.[9] Schwartz IB, Smith HL. Infinite subharmonic bifurcation in an SEIR epidemicmodel. J Math Biol 1983; 18:233-253.[10] Kamo M, Sasaki A. The effect of cross-immunity and seasonal forcing in amulti-strain epidemic model. Physica D 2002; 165:228-241.[11] London W, Yorke JA. Recurrent outbreaks of measles, chickenpox and mumps.I. Seasonal variation in contact rates. Am J Epidem 1973; 98(6):469-482.[12] Schenzle D. Age-structured model of pre- and post-vaccination measlestransmission. IMA J Math Appl Med Biol 1984; 1:169-191.[13] Buonomo B, Chitnis N, d’Onofrio A. Seasonality in epidemic models: aliterature review. Ricerche Mat 2018; 67:7–25.[14] London WP, Yorke JA. Recurrent outbreaks of measles, chikenpox and mumps.1. Seasonal variation in contact rates. Am J Epidemiol 1973; 98:453-468.
15] Fine PEM, Clarkson JA. Measles in England and Wales - I. An analysis offactors underlying seasonal patterns. Int J Epidemiol 1982; 11:5-14.[16] Keeling MJ, Grenfell BT. Desease extinction and community size: modeling thepersistence of measles. Science 1997; 275:65-67.[17] Earn DJ, Rohani P, Bolker BM, Grenfell BT. A simple model for complexdynamical transitions in epidemics. Science 2000; 287:667-670.[18] Stone L, Olinky R, Huppert A. Seasonal dynamics of recurrent epidemics.Nature 2007; 446:533-536.[19] Diedrichs DR, Isihara PA, Buursma DD. The schedule effect: can recurrentpeak infectious be reduced without vaccines, quarantines or school closings?Math Biosci 2014; 248:46-53.[20] Axelsen JB, Yaari R, Grenfell BT, Stone L. Multiannual forecasting of seasonalinfluenza dynamics reveals climatic and evolutionary drivers. Proc Natl AcadSci 2014; 111:9538-9542.[21] Billings L, Schwartz IB. Exiting chaos with noise: unexpected dynamics inepidemic outbreaks. J Math Biol 2002; 44:31-48.[22] Olsen LF, Schaffer WM. Chaos versus noisy periodicity: alternative hypothesesfor childhood epidemics. Science 1990; 249:499-504[23] Aron JL, Schwartz IB. Seasonality and period-doubling bifurcations in anepidemic model. J Theor Biol 1984; 110:665-679.[24] Keeling M, Rohani P, Grenfell BT. Seasonally forced disease dynamics exploredas switching between attractors. Physica D 2001; 148:317-335.[25] Yang J, Qu Z, Hu G. Duffing equation with two periodic forcings: the phaseeffect. Physical Review E 1996; 53:4402-4413.[26] Zambrano S, Allaria E, Brugioni S, Leyva I, Meucci R, Sanju´an MAF, ArecchiFT. Numerical and experimental exploration of phase control of chaos. Chaos2006; 16:1-9(013111).[27] Zambrano S, Seoane JM, Mari˜no IP,. Sanju´an MAF, Euzzor S, Meucci R,Arecchi FT. Phase control of excitable systems. New Journal of Physics 2008;10:1-12(073030).[28] Joseph SK, Mari˜no IP, Sanju´an MAF. Effect of the phase on the dynamics of aperturbed bouncing ball system. Commun Nonlinear Sci Numer Simulat 2012;17:3279-3286.[29] Used J, Wagemakers A, Sanju´an MAF. Regularization of map-based neuronmodels using phase control. Discontinuity, Nonlinearity and Complexity 2012;1:69-78.[30] Parker T, Chua LO. Practical numerical algorithms for chaotic systems.Springer-Verlag; 1989.
31] Greenman JV, Pasour VB. Threshold dynamics for periodically forcedecological systems: The control of population invasion and exclusion. Journalof Theoretical Biology 2012; 295:154-167.[32] Carlson CJ, Gomez ACR, Bansal S, Ryan SJ. Misconceptions about weatherand seasonality must not misguide COVID-19 response. Nat Commun 2020;11:4312.[33] Fan Wu, Su Zhao, Bin Yu, Yan-Mei Chen, Wen Wang, Zhi-Gang Song, Yi Hu,Zhao-Wu Tao, Jun-Hua Tian, Yuan-Yuan Pei, Ming-Li Yuan, Yu-Ling Zhang,Fa-Hui Dai, Yi Liu, Qi-Min Wang, Jiao-Jiao Zheng, Lin Xu, Edward C. Holmes,Yong-Zhen Zhang. A new coronavirus associated with human respiratory diseasein China. Nature 2020; 579:265-269.[34] Qianying Lin, Shi Zhao, Daozhou Gao, Yijun Lou, Shu Yang, Salihu S.Musa,Maggie H.Wang, Yongli Cai, Weiming Wang, Lin Yang, Daihai He. A conceptualmodel for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, Chinawith individual reaction and governmental action. International Journal ofInfectious Diseases 2020; 93:211-216.[35] Guan W et al. Early Transmission Dynamics in Wuhan, China, of NovelCoronavirus-Infected Pneumonia. N Engl J Med 2020; 382:1199-1207.[36] Zhu N et al. Clinical Characteristics of Coronavirus Disease 2019 in China, TheNew England Journal of Medicine 2020; 382:727.31] Greenman JV, Pasour VB. Threshold dynamics for periodically forcedecological systems: The control of population invasion and exclusion. Journalof Theoretical Biology 2012; 295:154-167.[32] Carlson CJ, Gomez ACR, Bansal S, Ryan SJ. Misconceptions about weatherand seasonality must not misguide COVID-19 response. Nat Commun 2020;11:4312.[33] Fan Wu, Su Zhao, Bin Yu, Yan-Mei Chen, Wen Wang, Zhi-Gang Song, Yi Hu,Zhao-Wu Tao, Jun-Hua Tian, Yuan-Yuan Pei, Ming-Li Yuan, Yu-Ling Zhang,Fa-Hui Dai, Yi Liu, Qi-Min Wang, Jiao-Jiao Zheng, Lin Xu, Edward C. Holmes,Yong-Zhen Zhang. A new coronavirus associated with human respiratory diseasein China. Nature 2020; 579:265-269.[34] Qianying Lin, Shi Zhao, Daozhou Gao, Yijun Lou, Shu Yang, Salihu S.Musa,Maggie H.Wang, Yongli Cai, Weiming Wang, Lin Yang, Daihai He. A conceptualmodel for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, Chinawith individual reaction and governmental action. International Journal ofInfectious Diseases 2020; 93:211-216.[35] Guan W et al. Early Transmission Dynamics in Wuhan, China, of NovelCoronavirus-Infected Pneumonia. N Engl J Med 2020; 382:1199-1207.[36] Zhu N et al. Clinical Characteristics of Coronavirus Disease 2019 in China, TheNew England Journal of Medicine 2020; 382:727.