# Seasonally Forced SIR Systems Applied to Respiratory Infectious Diseases, Bifurcations, and Chaos

Authors: Nico Stollenwerk; Stefano Spaziani; Javier Mar; Irati Eguiguren Arrizabalaga; Damián Knopoff; Nicole Cusimano; Vizda Anam; Akhil Shrivastava; Maíra Aguiar
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

Summary. We investigate models to describe respiratory diseases with fast mutating virus pathogens such that after some years the aquired resistance is lost and hosts can be infected with new variants of the pathogen. Such models were initially suggested for respiartory diseases like influenza, showing complex dynamics in reasonable parameter regions when comparing to historic empirical influenza like illness data, e.g., from Ille de France. The seasonal forcing typical for respiratory diseases gives rise to the different rich dynamical scenarios with even small parameter changes. Especially the seasonality of the infection leads for small values already to period doubling bifurcations into chaos, besides additional coexisting attractors. Such models could in the future also play a role in understanding the presently experienced COVID-19 pandemic, under emerging new variants and with only limited vaccine efficacies against newly upcoming variants. From first period doubling bifurcations, we can eventually infer at which close by parameter regions complex dynamics including deterministic chaos can arise.


---

## Body

## 1. Introduction

Some of the first successful modelling approaches in epidemiology which show complex dynamical behaviour were used to describe the observed large fluctuations in childhood diseases in historic prevaccination eras. Childhood diseases are typically highly contagious, but human hosts have developed life long immnunity against these pathogens, such that only children were typically infected but then protected for the remaining life time.But most endemic diseases are less infective, still showing at times complex dynamical patterns of numbers of infected in long-term time series recordings. Among those are respiratory diseases like influenza and other influenza-like illnesses (ILIs), where in winter typically larger or smaller outbreaks are observed over the years [1]. By only assuming a waning immunity of a few years and moderate infectivity, much below the ones observed for childhood diseases, we can describe already such large seasonal fluctuations with simple SIR-type models including a sinusoidal seasonal forcing [2].Here, we investigate such models for seasonal respiratory diseases with fast mutating pathogens in their capacity of generating rich bifurcation scenarios and deterministic chaos already for small variations of seasonality. Implications for the eventual future dynamical development of COVID-19 are outlined but will still depend on the escape capacity of new variants from already existing immunity against the initial variants either by natural infection or by vaccines developed with the knowledge of the initial variants.For data analysis and initial phase dynamics of COVID-19, especially in the Basque Country, seasonality was still difficult to quantify [3–5], but now after the first year and a half, it becomes clearer that there is a certain seasonal effect also in COVID-19, but now already overlayed by vaccination campaigns effects. Still, a lot of things changed in our knowledge about COVID-19, e.g., the intially high case fatality ratios reported in 2020 [6] are by now widely accepted to be lower than initially though, due to the wide discovery of asymptomatically or mildly infected [7]. Hence, some comparison to other respiratory diseases, including comparative data analysis, also might help here in the future evaluation of some of the key parameters such as the seasonality and the waning immunity/vaccine escape.


## 2. Seasonally Forced SIR Models

The seasonally forced SIR system can be given by the following reaction scheme:(1)S+I⟶βtI+I,I⟶γR,R⟶αS,S⟶ρI,and its corresponding dynamics of probabilities is discussed below. The population of size N is divided into susceptible individuals S, infected I, and recovered R, which can become susceptible again after a waning immunity period, given by the transition probability α. Further, β is the infection rate and γ the recovery rate of the disease. Further, ρ is an additional import, which leads susceptibles from the study population to become infected due to outside contacts. The seasonal forcing of the infection rate is given by
(2)βt=β0·1+θ·cosωt+φ,modelling the increased infectivity probability in colder weather in winter. Here, we first consider the deterministic skeleton, given by the mean field approximation of the underlying stochastic process [8]. The SIR model in its mean field ODE system, first without seasonal forcing but already including import ρ, is given by
(3)dSdt=αR−βNS·I+ρN,dIdt=βNS·I+ρN−γI,dRdt=γI−αR.Hence, we have a stationary state, without seasonality, given by
(4)R∗=γαI∗,(5)S∗=γβNI∗I∗+ρN,both still depending on I∗ and finally,
(6)I∗N=12αγ+α1−γβ−ρ+14αγ+α1−γβ−ρ2+αγ+α·ρ.This stationary state is a fixed point in the state space of the variablesS and I.Adding the seasonal forcing βt=β0·1+θ·cosωt+φ gives for very small seasonality θ a limit cycle around the fixed point and increases the state space by one more variable to a three-dimensional space of St, It, and βt. As parameter values, we consider influenza-like dynamical behaviour, for which due to mutuating influenza viruses the immunity gets lost after some years, hence giving on average a waning immunity of 6 years. Relevant paremeters are α=1/6y, γ=1/3d=365/3y−1, β0=1.5·γ, and lnρ=−15 and varying θ. For an analysis based on the time available influenza data from France, we took θ=0.12 as analyzed in detail in an earlier publication [2]. Here, we consider from the nonforced system, i.e., θ=0, and no import a gradual increase of the seasonality, and observe already for very small seasonality bifurcations and complex dynamics, much below θ=0.12.


### 2.1. Seasonal Forcing Included in the SIR System Leads to Bifurcations into Chaos

From an initial limit cycle around the nonseasonal fixed point, we already observe for a seasonality of, e.g.,θ=0.038, a deformed limit cycle of period length of two years (not of one year, the focing period length), see Figure 1(a). By increasing the seasonality further, to, e.g., θ=0.043, we observe a period doubling, from the two-year limit cycle to a four-year limit cycle, see Figure 1(b). In the time series of infected, this corresponds to a roughly two-year cycle with a large outbreak in one year and a small in the next season, but after that the large outbreak is not exactly repeated and also the low outbreak in the fourth season is slightly altered from the second season, and only after that the outbreak sizes are repeated.Figure 1
State space plots of period doubling. In state space projection of the variablesSt and It, with complete state space also including βt, we plot the limit cycles observed after cutting transient dynamics, for the following parameter values of the seasonality θ: (a) θ=0.038, (b) θ=0.043, and (c) θ=0.0446.
(a)(b)(c)Then, already in a very close by seasonality ofθ=0.0446, we observe a further priod doubling to a period length of 8 years, see Figure 1(c). The four-year cycle is not any more exactly matched but a second slightly altered turn occurs. This period doubling continues and leads to the period doubling rout to chaos, at which a deterministically chaotic attractor appears, with sensitivity to initial conditions, as, e.g., measured via positive Lyapunov exponents, see further below on the terms of deterministic chaos and Lyapunov exponents.Since we have initially a shift from a one-year limit cycle to a two-year limit cycle, the best to obtain good bifurcation diagrams is to consider stroboscopic maps, i.e., the value at each year at the same time of the year varied continuously with seasonalityθ, see Figure 2, where for low seasonality θ=0.038 two points appear, the value of low outbreak in one year and the value of high outbreak in the next. Then, we observe the period doubling, with four values for a fixed seasonality value, e.g., of θ=0.043, the value which corresponds to the state space plot in Figure 1(b).Figure 2
Bifurcation diagram around the period doubling. Now, stroboscopic values ofItn are used, e.g., each first of January. tn+1=tn+T with a period length of the seasonal forcing of T=1y. For graphical visibility, we plot the logarithm of the infected.Then, the sequential period doublings appear in the bifurcation diagram, Figure2, from which we can roughly read off the parameter values at which the bifurcations from one periodic limit cycle to the next double periodic limit cycle appear. We have as a first attempt the following values of θ2=0.0402, θ3=0.0441, and θ4=0.0449. From these first values, we can calculate a “speed of bifurcation accumulation” δ4=Δθ3/Δθ4=θ3−θ2/θ4−θ3=4.875. In practice, e.g., in the case of COVID-19, we could have a lower transmission in summer, and in the first few winter seasons, because of still occuring control measures due to mask wearing in closed space etc., we would have a dampened low seasonality. Then once such control measures will be relaxed, we could observe an increase in seasonality of infections, up to a point where we would observe a low ourbreak in one season and a higher outbreak in the next, burning out the remaining susceptibles, and a low season again, until the waning immunity gives rise to newly susceptibles via variant mutations of the pathogen. Ultimately, we could observe larger fluctuations as are, e.g., observed in seasonal influenza, for which we have a decade long outbreak data. Though, in general influenza, viruses are considered to mutate easier than the circulating human coronaviruses pre-COVID-19, the presently upcoming new variants leave some doubts about much milder circumstances in the case of the present COVID-19 epidemic to become also endemic with smaller and larger outbreak in future seasons.We will now apply some dynamical systems tools in order to understand better the implications of such “speed of bifurcation accumulation”, and especially the observation ofδ≈4.8 from our rough first attempt of quantification given above.


## 2.1. Seasonal Forcing Included in the SIR System Leads to Bifurcations into Chaos

From an initial limit cycle around the nonseasonal fixed point, we already observe for a seasonality of, e.g.,θ=0.038, a deformed limit cycle of period length of two years (not of one year, the focing period length), see Figure 1(a). By increasing the seasonality further, to, e.g., θ=0.043, we observe a period doubling, from the two-year limit cycle to a four-year limit cycle, see Figure 1(b). In the time series of infected, this corresponds to a roughly two-year cycle with a large outbreak in one year and a small in the next season, but after that the large outbreak is not exactly repeated and also the low outbreak in the fourth season is slightly altered from the second season, and only after that the outbreak sizes are repeated.Figure 1
State space plots of period doubling. In state space projection of the variablesSt and It, with complete state space also including βt, we plot the limit cycles observed after cutting transient dynamics, for the following parameter values of the seasonality θ: (a) θ=0.038, (b) θ=0.043, and (c) θ=0.0446.
(a)(b)(c)Then, already in a very close by seasonality ofθ=0.0446, we observe a further priod doubling to a period length of 8 years, see Figure 1(c). The four-year cycle is not any more exactly matched but a second slightly altered turn occurs. This period doubling continues and leads to the period doubling rout to chaos, at which a deterministically chaotic attractor appears, with sensitivity to initial conditions, as, e.g., measured via positive Lyapunov exponents, see further below on the terms of deterministic chaos and Lyapunov exponents.Since we have initially a shift from a one-year limit cycle to a two-year limit cycle, the best to obtain good bifurcation diagrams is to consider stroboscopic maps, i.e., the value at each year at the same time of the year varied continuously with seasonalityθ, see Figure 2, where for low seasonality θ=0.038 two points appear, the value of low outbreak in one year and the value of high outbreak in the next. Then, we observe the period doubling, with four values for a fixed seasonality value, e.g., of θ=0.043, the value which corresponds to the state space plot in Figure 1(b).Figure 2
Bifurcation diagram around the period doubling. Now, stroboscopic values ofItn are used, e.g., each first of January. tn+1=tn+T with a period length of the seasonal forcing of T=1y. For graphical visibility, we plot the logarithm of the infected.Then, the sequential period doublings appear in the bifurcation diagram, Figure2, from which we can roughly read off the parameter values at which the bifurcations from one periodic limit cycle to the next double periodic limit cycle appear. We have as a first attempt the following values of θ2=0.0402, θ3=0.0441, and θ4=0.0449. From these first values, we can calculate a “speed of bifurcation accumulation” δ4=Δθ3/Δθ4=θ3−θ2/θ4−θ3=4.875. In practice, e.g., in the case of COVID-19, we could have a lower transmission in summer, and in the first few winter seasons, because of still occuring control measures due to mask wearing in closed space etc., we would have a dampened low seasonality. Then once such control measures will be relaxed, we could observe an increase in seasonality of infections, up to a point where we would observe a low ourbreak in one season and a higher outbreak in the next, burning out the remaining susceptibles, and a low season again, until the waning immunity gives rise to newly susceptibles via variant mutations of the pathogen. Ultimately, we could observe larger fluctuations as are, e.g., observed in seasonal influenza, for which we have a decade long outbreak data. Though, in general influenza, viruses are considered to mutate easier than the circulating human coronaviruses pre-COVID-19, the presently upcoming new variants leave some doubts about much milder circumstances in the case of the present COVID-19 epidemic to become also endemic with smaller and larger outbreak in future seasons.We will now apply some dynamical systems tools in order to understand better the implications of such “speed of bifurcation accumulation”, and especially the observation ofδ≈4.8 from our rough first attempt of quantification given above.


## 3. Seasonally Forced ODE Models to Maxima Return Maps and Its Simplifications

When considering a seasonally forced SIR system in its chaotic regime, in any parameter combination as long as it is chaotic, we can take a “maxima return map”, i.e., the maximum at each outbreak and plot this against the next maximum. Alternatively, Poincaré sections in the state space could be considered, or the abovementioned stroboscopic maps. A typical maxima return map is given in Figure3(a). After a very large outbreak, we observe rather small outbreaks next, a large In followed by a small In+1. For intermediate small In, we observe rather large In+1. But for very small In, also, the In+1 is rather small, since the susceptibles are not yet build up sufficiently to generate a next large outbreak. Finally, intermediate outbreaks are often followed by another intermediate outbreak.Figure 3
(a) Maxima return map of a seasonally forced SIR model in its chaotic parameter region. (b) Simplified one-dimesional version of the maxima retun map, the Ricker map, here in its simples functional form, not adjusted for numerical comparison with the maxima return map.
(a)(b)Though such discrete maps generated from 3-dimensional ODE flows have to be two dimensional, as seen in Figure3(a), the described sequence of In+1 following In can be roughly represented by a simple functional form of x·e−x. Hence, we can give a well know map from ecology, the Ricker map, in its implest form
(7)xn+1=fxn=axne−xn,as a crude description of the dynamical behaviour of the maxima return map originated from a seasonally forced SIR system. Though the maxima return map was obtained in the chaotic region of the seasonal SIR model, also in the sketchy Ricker map, we observe period doubling bifurcations for certain parameter values a, see Figure 4(a), giving the bifurcation diagram of the Ricker map for, e.g., a∈ [4, 9].Figure 4
(a) Bifurcation diagram for the Ricker map. (b) Lyapunov spectrum for the Ricker map in the same parameter region for comparison.
(a)(b)The period doubling bifurcations cumulate into a fuzzy clowd of points in the bifurcation diagram, Figure4(a). This clowd is the onset of deterministically chaotic attractors. Such deterministic chaos can be quantified by measuring the exponential convergence or divergence along any trajectory in its time average after very long time, the so called “Lyapunov exponent”
(8)λ=limN⟶∞1N∑n=1Nlndfdxxn,which can be plotted with varying parameter a to give a Lyapunov spectrum, see Figure 4(b). We use the same range of the parameter a in Figure 4(b) as in Figure 4(a) to be able to compare the bifurcation diagram and Lyapunov spectrum.The limit cycles have typically exponential convergence, hence negative Lyapunov exponents. Only at the bifurcation points at each period doubling, one periodic cycle becomes unstable, and the double periodic cycle becomes stable. Hence, stability changes at the bifurcation points, and the Lyapunov exponents are zero there, going negative again for the next double limit cycle. Such characterization of bifurcations via zero Lyapunov exponents could, e.g., be used in graphical two-dimensional bifurcation diagrams [10].After many observed, actually infinitely many, period doubling bifurcations, we observe positive Lyapunov exponents, from arounda≈14.2 on, indicating that the period doubling cascade culminates in chaotic attractors. Again, we can measure the speed of period doubling bifurcations into a new dynamic regime, the chaotic attractors, and find numerically a similar value as our rough first attempt in the seasonal SIR system gave.Namely, we determine subsequent bifurcation points from zooming into the bifurcation diagram in smaller and smaller regions of the parametera, see, e.g., Figure 5 for the final step, and observe ac1=e2=7.389056099, which can be calculated analytically, then ac2=12.509, giving a difference of Δa2≔ac2−ac1=5.12.Figure 5
Last step zooming into the bifurcation diagram to determinac6=14.76137, where the value is read off from the underlying data file.Continuing this analysis further, we obtainac3=14.2443, Δa3≔ac3−ac2=1.735, then ac4=14.6527, Δa4=0.4084, and ac5=14.7421, Δa5=0.0894, and finally ac6≈14.76137, Δa6=0.0193. Hence, for the speed of bifurcations, we obtain δ3≔Δa2/Δa3=5.12/1.735=2.951, δ4≔Δa3/Δa4=1.735/0.4084=4.248, δ5≔Δa4/Δa5=0.4084/0.0894=4.5682, and δ6≔Δa5/Δa6=0.0.0894/0.0193=4.6321. Hence, we move towards a δ∞=4.6692, which gives for subsequent bifurcation points ac,j+1≈ac,j+1/δ·ac,j−ac,j−1, ones δj+1≈δj, measuring the convergence speed into chaos at ac,∞≈14.2.Remember that we now found values of the speed of bifurcations close toδ∞=4.6692 in the seasonal SIR model as well as in the much simplified Ricker map, via the maxima return map. This already indicates a certain “universality” of δ∞=4.6692. If we would refine our analysis of the SIR-system in the same way we did for the Ricker map, we would expect similar values of δj towards δ∞=4.6692.The numerical agreement, here, demonstrated only roughly, can be improved and holds due to the so called “Feigenbaum conjecture,” namely, that the speed in period doubling bifurcation sequences into chaos ultimately, close to chaos, the numerically same value is obtained equally in ODE systems and in maps, one-dimensional or higher dimensional. A proof of the Feigenbaum constant is typically obtained for one-dimensional maps on the unity interval, as sketched below, but holds empirically much wider, for ODE systems and experimental physical systems [11], and here also shown for epidemiological models.


## 4. Feigenbaum Renomalization Gives Analytic Handle on Period Doubling Convergence Speed

We have obtained approximations ofδ∞=4.6692 from the bifurcation diagrams of different systems. The fastest way to obtain good approximations of δ∞ is via the so-called Feigenbaum renormalization [12], which resembles renormalizations in statistical systems close to criticality [8, 13] and renormalization in quantum physics [13]. Originally, Feigenbaum used a simple quadratic map on the unit interval, namely,
(9)xn+1=fxn=axn1−xn,and derived from this a functional equation of exact selfsimilarity at the renormalization fixed point [12], the so-called Cvitanović equation, which however is most widely known from a linear transformed version of the quadratic map, now on the interval between -1 and +1. The solution of this Cvitanović equation and its linearization around the fixed point [14] give the numerically best values of δ∞≕δF, the universal Feigenbaum constant δF, in which the value was however already known before the publications of Feigenbaum [15].To describe the self-similarity in the Feigenbaum map, we first iterate the map for a moderate parameter value, e.g.,a=2.9, and observe that the trajectory xn converges to a fixed point, but for a value of, e.g., a=3.3 keeps oscillating in a period 2, see Figure 6(e).Figure 6
Self-similarity of the Feigenbaum iteration. Period 8 limit cycle fora=3.554 in (a) and in fourth iteration giving a period 2 limit cycle in (b). Period 4 for a=3.5 in (c) and second iterate in (d). Period 2 for a=3.3 in (e) and self-similarity of this in (f) zooming into the previous higher iterate in (b) and inverted and zooming in (d).
(a)(b)(c)(d)(e)(f)Then, for, e.g.,a=3.5, we find a period 4 cycle in the Feigenbaum map, Figure 6(c). If we iterate the map, i.e., xn+2=ffxn≕f2xn, the period 4 becomes a period 2 limit cycle again, see Figure 6(d). And for, e.g., a=3.554, we find a period 8 limit cycle in the original Feigenbaum map, Figure 6(a), which becomes a period 2 limit cycle in the fourth iterate xn+4=f4xn, Figure 6(b). Now, if we compare the fourth iterate in Figure 6(b) with the first iterate in Figure 6(f), we find the same graphical features from the first iterate again in the fourth iterate by simply zooming into the center of the unit interval, around xn=0.5, compare Figures 6(f) and 6(b).Now, the original map Figure6(f) can be also compared with the second iterate Figure 6(d), when we flip the inner part around xn=0.5 and zoom in or mathematically rescale x-axis and y-axis. This rescaling can be given mathematically in the following way:
(10)1−fai+1fai+11−x≈1α·faiα·x,with a scaling factor αF (or here simply α since it cannot be confused with the parameter α in the SIR models) and flipping via 1−x in the x-axis and 1−f2 in the y-axis. Though this approximate self-similarity would be sufficient for the subsequent analysis, a map just being linearly transformed into the interval -1 to +1, namely, xn+1=1−μxn2≕fμxn is generally used with the approximate self-similarity
(11)−αfμi+1fμi+1−xα≈fμix,with scale factor again called α. This defines an operator T which maps fμix into fμi+1(12)Tfx≔−αff−xα,with the generally assumed assumption in renormalization, that after infinitely many iterations, we arrive at the critical point of the onset of chaos to an exact self-similarity; hence, T has a fixed point function gx for period doubling going to infinity fulfilling
(13)Tgx=gx,which gives explicitly the so-called Cvitanović equation
(14)gx=−αgg−xαfor the universal constant α=2.5029 and universal function gx. The constant α=2.5029 is the first Feigenbaum constant, while however we are here more interested in the second Feigenbaum constant δF. The renormalization theory now gives an attractive high-dimensional manifold, attracting towards the fixed point Tgx=gx, and one unstable direction in function space with dominating eigenvalue, which turns out to be δF we are looking for [16].Hence, we investigate the linearization ofT around the renormalization fixed point Tg=g by
(15)Lhx≔δTfxδfxgx=Tgx+hx−Tgx,for small perturbation hx, resulting explicitly in
(16)Lhx=−αg′g−xα·h−xα+hg−xα,and obtain the eigenvalue/eigenfunction problem in function space
(17)Lhx=δF·hx,with identified
(18)δF≔limi⟶∞ai−ai−1ai+1−ai=4.6692,at the onset of chaos, after infinitely many period doubling bifurcations, at parameter value a with value a∞=3.5699456.The fixed point equationTgx=gx then can be solved via popynomial ansatz of gx giving also the value for α and the linearization Lhx=δF·hx giving in addition δF, each in respective accuracy depending on the length of the considered polynomials. We will only give a brief glance at the first steps below, which however already give relatively good approximations for the two Feigenbaum constants α and δF, with further refinements easily obtained to high precision [17].


### 4.1. Practical Solution of the Functional Equations

We start with the polynomial Ansatz for the universal functiong quadratically via gx=g0+g1x+g2x2+Ox3 to be inserted into Tgx=gx=−αgg−x/α. This gives via coefficient comparison in orders of x the results g1=0, g0=1 and g2=−α/2, and from this
(19)α≈1+3=2.73,compared to the higher order polynomial results of
(20)α=2.5029,eventually with much higher precision [17] than reported here. It is however surprising that already the here discussed approximation give the correct order of magnitude of α.Now,hx is inserted into Lhx=δF·hx which already gives in zeroth order
(21)δF≈α2−α=3+3=4.73,compared to higher order result of
(22)δF=4.6692,again in good agreement in orders of magnitude. And of course from here, we can continue in higher orders of polynomials further.Though in the renormalization flow the diverging function space direction is given byδF, in terms of bifurcation accumulation δF measures how quickly complexity increases in terms of bifurcations with higher and higher orders and soon reaches deterministically chaotic attractors, due to the relatively high value of δF. Further, α measures the amplitude increase between bifurcations, here in the epidemiological examples the maximal numbers of infected in each year’s outbreaks.Not only that, the Feigenbaum constantsα and δ are found empirically in many period doubling routes to chaos in ODE systems and in physical systems (laser experiments, etc.) [11], hence universal constants, we also observe them in epidemiological systems, as indicated here. Hence, one could observe the first indications of increasing seasonality giving rise to increased complexity in time series, i.e., from periodic oscillations to double periodic oscillations, one would be in the position to expect more complex dynamics with smaller and smaller increases of seasonality. We will now give a first look into stochastic versions of the epidemiological models already in the chaotic parameter region, still observing that to some extend the notions of the underlying deterministic skeleton help in understanding the system’s behavious, though enhancing the fluctuations by noise.


## 4.1. Practical Solution of the Functional Equations

We start with the polynomial Ansatz for the universal functiong quadratically via gx=g0+g1x+g2x2+Ox3 to be inserted into Tgx=gx=−αgg−x/α. This gives via coefficient comparison in orders of x the results g1=0, g0=1 and g2=−α/2, and from this
(19)α≈1+3=2.73,compared to the higher order polynomial results of
(20)α=2.5029,eventually with much higher precision [17] than reported here. It is however surprising that already the here discussed approximation give the correct order of magnitude of α.Now,hx is inserted into Lhx=δF·hx which already gives in zeroth order
(21)δF≈α2−α=3+3=4.73,compared to higher order result of
(22)δF=4.6692,again in good agreement in orders of magnitude. And of course from here, we can continue in higher orders of polynomials further.Though in the renormalization flow the diverging function space direction is given byδF, in terms of bifurcation accumulation δF measures how quickly complexity increases in terms of bifurcations with higher and higher orders and soon reaches deterministically chaotic attractors, due to the relatively high value of δF. Further, α measures the amplitude increase between bifurcations, here in the epidemiological examples the maximal numbers of infected in each year’s outbreaks.Not only that, the Feigenbaum constantsα and δ are found empirically in many period doubling routes to chaos in ODE systems and in physical systems (laser experiments, etc.) [11], hence universal constants, we also observe them in epidemiological systems, as indicated here. Hence, one could observe the first indications of increasing seasonality giving rise to increased complexity in time series, i.e., from periodic oscillations to double periodic oscillations, one would be in the position to expect more complex dynamics with smaller and smaller increases of seasonality. We will now give a first look into stochastic versions of the epidemiological models already in the chaotic parameter region, still observing that to some extend the notions of the underlying deterministic skeleton help in understanding the system’s behavious, though enhancing the fluctuations by noise.


## 5. Stochastic SIR Systems

After demonstrating that there is a route to chaos in such simple seasonal SIR systems describing respiratory diseases, we now investigate the robustness against stochastic influences, in a case study already well inside the parameter region of expected deterministically chaotic dynamics. Sometimes, there are still concerns expressed that complex dynamics rather might be decreased by noise, but on the contrary, we observe that taking noise appropriately into account including state dependence, we observe even under large population sizes complex behaviour as to be expected from the previous deterministic process analysis of mean field dynamics.


### 5.1. Stochastic Versions of the Seasonally Forced SIR Model

From reaction scheme 1, we can immediately give the master equation for the SIR system as dynamics of the probabilities(23)ddtpS,I,t=βtNS+1I−1+ρNpS+1,I−1,t+γI+1pS,I+1,t+αN−S−1−IpS−1,I,t−βtNSI+ρN+γI+αN−S−IpS,I,t.The master equation can be written in a generic form using densitiesx1≔S/N and x2≔I/N instead of the total numbers of individuals in the population classes S and I, hence state vector x¯≔x1,x2tr, as
(24)ddtpx¯,t=∑j=1nNwjx¯+Δx¯j⋅px¯+Δx¯j,t−Nwjx¯⋅px¯,t,with n=3 different transitions and small deviation from state x¯ as Δx¯j≔1/N·r¯j.For the master equation in densitiesx1≔S/N and x2≔I/N, we have now transitions wjx¯ and its vectors r¯j given by
(25)w1x¯=βtx1x2+ρ,r¯1=1,−1tr,w2x¯=γx2,r¯2=0,1tr,w3x¯=α1−x1−x2,r¯3=−1,0tr..We can perform simulations of the master equation directly using the Gillespie algorithm [18, 19], observing that for large population sizes N the simulations become very slow. So, we will use the Kramers-Moyal expansion to obtain a Fokker-Planck equation as approximation of the master equation, for which the corresponding stochastic differential equation system can be simulated much faster [9, 20, 21].


### 5.2. Fokker-Planck Approximation of the SIR Model

From the master equation in densities, Equation (24) for the SIR system, we use Taylor’s expansion
(26)wjx¯+Δx¯j·px¯+Δx¯j,t=∑ν=0∞1ν!Δx¯j·∇x¯νwjx¯px¯,t,giving to second order in 1/N a Fokker-Planck equation
(27)∂∂tpx¯,t=−∇x¯∑j=1n−r¯j⋅wjx¯px¯,t+σ22∑j=1nr¯j⋅∇x¯2wjx¯px¯,t,with
(28)∇x¯=∂∂x1∂∂x2=∂x¯,or in a different notation
(29)∂∂tpx¯,t=−∂x¯f¯x¯px¯,t+σ22∂⟶x¯G2x¯,t∂⟵x¯,using simply a quadratic form ∂x¯⟶G2x¯px¯,t∂x¯⟵ here with
(30)∂x¯⟶G2p∂x¯⟵=∂∂x1,∂∂x2·g11g12g21g222px¯,t·∂∂x1⟵∂∂x2⟵,(31)f¯x¯=∑j=1nf¯jx¯=∑j=1n−r¯j·wjx¯,G2x¯=∑j=1nGj2x¯=∑j=1nr¯j·r¯jtrwjx¯.The Fokker-Planck equation gives a stochastic differential equation system withσ=1/N, and in the SIR model, the two-dimensional Gaussian normal noise vector ε¯t=εx1t,εx2ttr as
(32)ddtx¯=f¯x¯+σGx¯·ε¯t,and using matrix square root from eigenvalue-eigenvector decomposition G2x¯=TΛT−1 as Gx¯=TΛTtr to be numerically implemented easily, and much faster than the Gillespie algorithm for the master equation [18, 19]. To further speed up the SDE simulations, we can relax the condition of G being a quadratic matrix and allow a matrix B with BBtr=G2, such that the expressions in B can become quite simple and fast to compute [22]. Explicitly, we have for B the expression in the SIR model given as B=w1r¯1,w2r¯2,w3r¯3,, giving again G2=BBtr=∑j=1nr¯j·r¯jtrwjx¯ as required.In Figure7, we show simulations for the seasonal SIR model with the previously stated parameters and high seasonality θ=0.12, for the mean field ODE system, the master equation simulated via the Gillespie algorithm, and the stochastic differential equation systems using the quadratic G matrix and the nonquadratic B matrix for the covariances. We start in all four cases with the same initial conditions, hence, can observe the only slowly increasing differences of the stochastic processes from the mean field ODE system, in Figure 7(a) for the numbers of infected and in (b) for the susceptibles. In (c), we show in state space projection the deterministically chaotic attractor of the mean field ODE system, and in (d) also, the comparison with the stochastic process realizations, which preserve the state space structure, but eventually explore larger excursions into state space regions with small contraction to the deterministic attractor.Figure 7
Stochastic fluctuations around the chaotic attractor. Time series and state space plots of mean field ODE system, Gillespie algorithmic realization of master equation, and Fokker-Planck bases SDE realizations. In yellow, the deterministic skeleton of the mean field solution, see, e.g., in (c) its state space behaviour, and larger fluctuations in the stochastic realizations in (d). Note that we simulate the equations with high time resolution, but then sample with lower resolution, according to the schemes one would record empirical data, so that we see piecewise linear time evolutions graphically.
(a)(b)(c)(d)While we previously investigated the seasonal SIR system in this parameter region further in respect to parameter estimation [2] via the so-called η-ball method [23], we had besides the Gillespie algorithm only the multinomial approximation of the stochastic process at hand, both becoming equally slow for large population sizes [2], but now newly tested faster procedures based on the Fokker-Planck approximation to stochastic differential equation systems are available which work well for large population sizes outside the extinction boundary, as more recently explored in other epidemiological systems [24], including model comparison on the basis of Bayes factors [25].In the chemical literature, further aspects of speeding up the simulations from the classical Gillespie algorithm of exponential waiting times to the next reaction, also called Gillespie’s direct method, have been explored. e.g., close to the exponential waiting time algorithm, one can approximate in small time intervalsΔt the number of reactions by Poisson distributions, assuming little changes in the transitions wx¯ during the interval Δt, the so called τ-leap algorithms, which tend for increased Δt from the Poisson distributions again to the Gaussian distributions of the stochastic differential equation systems, as given in Equation (32) with the B-matrix approach, see, e.g., a detailed discussion by Gillespie [26], where the SDEs are referred to as chemical Langevin equations (CLE), parallel to the classical Langevin equation known in physics [9, 20]. Between exponential waiting times and τ-leap, another method, balancing speed of simulation with accuracy in stochastic trajectories still close to the exact direct method, is the so called kα method to evaluate waiting times after the kth event of some characteristic reactions, eventually strengthening the wx¯≈const. condition [26]. Here, the Γ-distribution as continuous time limit of discrete time Erlang k-distributions comes to play.However, the Gaussian approximation to the SDE systems, or Langevin equations, is of main importance, as among others, Gillespie emphasizes in 2001 as “The Langevin method therefore plays the important conceptual role of showing the stochastic simulation methods are related to the deterministic reaction rate equations RRE of traditional chemical kinetics.” [26]. This links the stochastic processes with our considerations of the ODE systems analyzed above, Section 2 and further.Other schemes combining the above with ideas of evaluating the state-dependent stochastic part of the Langevin equation in the mid of the intervalΔt, which naturally lead to Stratonovich versions of SDEs rather than Ito versions as discussed here, but then corrected with the Ito drift, namely, schemes of Milshtein type [9] are considered more recently [27], allowing larger time intervals Δt during the simulations. Basically, in the simplest case of a one-dimensional Langevin equation of the typ dx/dt=fx+σgxεt, the state-dependent noise term gx is evaluated in the interval xt,xt+Δt as gxt+c·xt+Δt−xt with c∈0,1, giving in its simplest version for c=0 the Ito version and for c=1/2 the Startonovich version [28], and adding to the drift fx the correction to Ito type, namely f−σ2/2g′g following obviously from the Fokker-Planck formulation [9, 20], gives the Milshtein method, either kept implicit to desired order in Δx≔xt+Δt−xt or made explicit like in the original Stratonovich version [28–30]. Such fully implicit stochastic methods [27] seem promising for highly accurate stochastic paths of chemical systems, speeding up from the original Gillespie direct method. Since we consider here applications to chaotic attractors, hence, the dynamics in state space being overall contracting, we are interested in typical trajectories close to the deterministic counterparts, and not yet consider crossing of boundaries between basins of attraction of coexisting attractors, as, e.g., analyzed in maps, derived from forced SIR-type models in other applications [31], so that higher precision trajectories’ simulations only would be beneficial, ones the basic concepts and future applications to respiratory disease dynamics will be settled.The present study shows that noise does not inhibit the complexity of the original determistic ODE systems, or in chemical literature terms, the reaction rate equations, but rather presevers or eventually might further enhance the complexity of the systems under study. Since we donot know yet in which parameter region a future COVID-19 pandemic will evolve with new variants of concern and reduction of vaccine efficacies, which leads to waning immunity, and on top, the exact seasonality of COVID-19 is still under large debate with rising numbers even in autumn and winter 2021 in already relatively well vaccinated populations, more detailed analyses of such systems of the type of seasonal SIR models have to rely on previous experiences in other respiratory diseases, which already show large fluctuations [32, 33]. Finally, empirical epidemiological data will decide to which extension refined numerical schemes will be needed, and on the other hand, bifurcation analyses and understanding of qualitative state space behaviours are needed to understand the systems under investigation [31].


## 5.1. Stochastic Versions of the Seasonally Forced SIR Model

From reaction scheme 1, we can immediately give the master equation for the SIR system as dynamics of the probabilities(23)ddtpS,I,t=βtNS+1I−1+ρNpS+1,I−1,t+γI+1pS,I+1,t+αN−S−1−IpS−1,I,t−βtNSI+ρN+γI+αN−S−IpS,I,t.The master equation can be written in a generic form using densitiesx1≔S/N and x2≔I/N instead of the total numbers of individuals in the population classes S and I, hence state vector x¯≔x1,x2tr, as
(24)ddtpx¯,t=∑j=1nNwjx¯+Δx¯j⋅px¯+Δx¯j,t−Nwjx¯⋅px¯,t,with n=3 different transitions and small deviation from state x¯ as Δx¯j≔1/N·r¯j.For the master equation in densitiesx1≔S/N and x2≔I/N, we have now transitions wjx¯ and its vectors r¯j given by
(25)w1x¯=βtx1x2+ρ,r¯1=1,−1tr,w2x¯=γx2,r¯2=0,1tr,w3x¯=α1−x1−x2,r¯3=−1,0tr..We can perform simulations of the master equation directly using the Gillespie algorithm [18, 19], observing that for large population sizes N the simulations become very slow. So, we will use the Kramers-Moyal expansion to obtain a Fokker-Planck equation as approximation of the master equation, for which the corresponding stochastic differential equation system can be simulated much faster [9, 20, 21].


## 5.2. Fokker-Planck Approximation of the SIR Model

From the master equation in densities, Equation (24) for the SIR system, we use Taylor’s expansion
(26)wjx¯+Δx¯j·px¯+Δx¯j,t=∑ν=0∞1ν!Δx¯j·∇x¯νwjx¯px¯,t,giving to second order in 1/N a Fokker-Planck equation
(27)∂∂tpx¯,t=−∇x¯∑j=1n−r¯j⋅wjx¯px¯,t+σ22∑j=1nr¯j⋅∇x¯2wjx¯px¯,t,with
(28)∇x¯=∂∂x1∂∂x2=∂x¯,or in a different notation
(29)∂∂tpx¯,t=−∂x¯f¯x¯px¯,t+σ22∂⟶x¯G2x¯,t∂⟵x¯,using simply a quadratic form ∂x¯⟶G2x¯px¯,t∂x¯⟵ here with
(30)∂x¯⟶G2p∂x¯⟵=∂∂x1,∂∂x2·g11g12g21g222px¯,t·∂∂x1⟵∂∂x2⟵,(31)f¯x¯=∑j=1nf¯jx¯=∑j=1n−r¯j·wjx¯,G2x¯=∑j=1nGj2x¯=∑j=1nr¯j·r¯jtrwjx¯.The Fokker-Planck equation gives a stochastic differential equation system withσ=1/N, and in the SIR model, the two-dimensional Gaussian normal noise vector ε¯t=εx1t,εx2ttr as
(32)ddtx¯=f¯x¯+σGx¯·ε¯t,and using matrix square root from eigenvalue-eigenvector decomposition G2x¯=TΛT−1 as Gx¯=TΛTtr to be numerically implemented easily, and much faster than the Gillespie algorithm for the master equation [18, 19]. To further speed up the SDE simulations, we can relax the condition of G being a quadratic matrix and allow a matrix B with BBtr=G2, such that the expressions in B can become quite simple and fast to compute [22]. Explicitly, we have for B the expression in the SIR model given as B=w1r¯1,w2r¯2,w3r¯3,, giving again G2=BBtr=∑j=1nr¯j·r¯jtrwjx¯ as required.In Figure7, we show simulations for the seasonal SIR model with the previously stated parameters and high seasonality θ=0.12, for the mean field ODE system, the master equation simulated via the Gillespie algorithm, and the stochastic differential equation systems using the quadratic G matrix and the nonquadratic B matrix for the covariances. We start in all four cases with the same initial conditions, hence, can observe the only slowly increasing differences of the stochastic processes from the mean field ODE system, in Figure 7(a) for the numbers of infected and in (b) for the susceptibles. In (c), we show in state space projection the deterministically chaotic attractor of the mean field ODE system, and in (d) also, the comparison with the stochastic process realizations, which preserve the state space structure, but eventually explore larger excursions into state space regions with small contraction to the deterministic attractor.Figure 7
Stochastic fluctuations around the chaotic attractor. Time series and state space plots of mean field ODE system, Gillespie algorithmic realization of master equation, and Fokker-Planck bases SDE realizations. In yellow, the deterministic skeleton of the mean field solution, see, e.g., in (c) its state space behaviour, and larger fluctuations in the stochastic realizations in (d). Note that we simulate the equations with high time resolution, but then sample with lower resolution, according to the schemes one would record empirical data, so that we see piecewise linear time evolutions graphically.
(a)(b)(c)(d)While we previously investigated the seasonal SIR system in this parameter region further in respect to parameter estimation [2] via the so-called η-ball method [23], we had besides the Gillespie algorithm only the multinomial approximation of the stochastic process at hand, both becoming equally slow for large population sizes [2], but now newly tested faster procedures based on the Fokker-Planck approximation to stochastic differential equation systems are available which work well for large population sizes outside the extinction boundary, as more recently explored in other epidemiological systems [24], including model comparison on the basis of Bayes factors [25].In the chemical literature, further aspects of speeding up the simulations from the classical Gillespie algorithm of exponential waiting times to the next reaction, also called Gillespie’s direct method, have been explored. e.g., close to the exponential waiting time algorithm, one can approximate in small time intervalsΔt the number of reactions by Poisson distributions, assuming little changes in the transitions wx¯ during the interval Δt, the so called τ-leap algorithms, which tend for increased Δt from the Poisson distributions again to the Gaussian distributions of the stochastic differential equation systems, as given in Equation (32) with the B-matrix approach, see, e.g., a detailed discussion by Gillespie [26], where the SDEs are referred to as chemical Langevin equations (CLE), parallel to the classical Langevin equation known in physics [9, 20]. Between exponential waiting times and τ-leap, another method, balancing speed of simulation with accuracy in stochastic trajectories still close to the exact direct method, is the so called kα method to evaluate waiting times after the kth event of some characteristic reactions, eventually strengthening the wx¯≈const. condition [26]. Here, the Γ-distribution as continuous time limit of discrete time Erlang k-distributions comes to play.However, the Gaussian approximation to the SDE systems, or Langevin equations, is of main importance, as among others, Gillespie emphasizes in 2001 as “The Langevin method therefore plays the important conceptual role of showing the stochastic simulation methods are related to the deterministic reaction rate equations RRE of traditional chemical kinetics.” [26]. This links the stochastic processes with our considerations of the ODE systems analyzed above, Section 2 and further.Other schemes combining the above with ideas of evaluating the state-dependent stochastic part of the Langevin equation in the mid of the intervalΔt, which naturally lead to Stratonovich versions of SDEs rather than Ito versions as discussed here, but then corrected with the Ito drift, namely, schemes of Milshtein type [9] are considered more recently [27], allowing larger time intervals Δt during the simulations. Basically, in the simplest case of a one-dimensional Langevin equation of the typ dx/dt=fx+σgxεt, the state-dependent noise term gx is evaluated in the interval xt,xt+Δt as gxt+c·xt+Δt−xt with c∈0,1, giving in its simplest version for c=0 the Ito version and for c=1/2 the Startonovich version [28], and adding to the drift fx the correction to Ito type, namely f−σ2/2g′g following obviously from the Fokker-Planck formulation [9, 20], gives the Milshtein method, either kept implicit to desired order in Δx≔xt+Δt−xt or made explicit like in the original Stratonovich version [28–30]. Such fully implicit stochastic methods [27] seem promising for highly accurate stochastic paths of chemical systems, speeding up from the original Gillespie direct method. Since we consider here applications to chaotic attractors, hence, the dynamics in state space being overall contracting, we are interested in typical trajectories close to the deterministic counterparts, and not yet consider crossing of boundaries between basins of attraction of coexisting attractors, as, e.g., analyzed in maps, derived from forced SIR-type models in other applications [31], so that higher precision trajectories’ simulations only would be beneficial, ones the basic concepts and future applications to respiratory disease dynamics will be settled.The present study shows that noise does not inhibit the complexity of the original determistic ODE systems, or in chemical literature terms, the reaction rate equations, but rather presevers or eventually might further enhance the complexity of the systems under study. Since we donot know yet in which parameter region a future COVID-19 pandemic will evolve with new variants of concern and reduction of vaccine efficacies, which leads to waning immunity, and on top, the exact seasonality of COVID-19 is still under large debate with rising numbers even in autumn and winter 2021 in already relatively well vaccinated populations, more detailed analyses of such systems of the type of seasonal SIR models have to rely on previous experiences in other respiratory diseases, which already show large fluctuations [32, 33]. Finally, empirical epidemiological data will decide to which extension refined numerical schemes will be needed, and on the other hand, bifurcation analyses and understanding of qualitative state space behaviours are needed to understand the systems under investigation [31].


## 6. Conclusions

We have investigated models able to describe respiratory disease outbreak dynamics, in cases relevant for influenza-like illnesses, and outlined implications for the future of the presently observed COVID-19 pandemic. Already very small seasonal forcing, typical for respiratory diseases, can generate complex dynamics including period doubling bifurcations into chaotic dynamics, which also shows to be robust under stochastic modelling versions. The Feigenbaum universality indicates very rapid increase of complexity with only small changes of seasonality.In initial studies, we also observed coexistences between different attractors, as also observable in seasonal SIR models in completely different parameter regions, describing childhood diseases [31], but the interplay between such co-existing attractors in the case of the present parameter regimes of respiratory diseases is work in progess and will be reported elsewhere.On the methodological side, results of renormalization theories can be applied in actual epidemiological systems, where previously critical fluctuations in COVID-19 dynamics were described [34], with universal critical exponents, which are analyzed in their universality by the renormalization theory in statistical physics [8, 13], a predecessor of the here described Feigenbaum renormalization [11]. Still such fruitfull techniques are not yet widely considered in biomathematics in their application aspects, but merely mentioned at times out of curiosity, e.g., the Feigenbaum map was described as a loose biological model for insect populations [35], respectively, the Ricker map [36] without consideration of its closeness to, e.g., maxima return maps of epidemiological systems.



---

*Source: 3556043-2022-03-03.xml*