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

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.


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 influenzalike 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][4][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.

Seasonally Forced SIR Models
The seasonally forced SIR system can be given by the following reaction scheme: 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 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 Hence, we have a stationary state, without seasonality, given by both still depending on I * and finally, This stationary state is a fixed point in the state space of the variables S 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 SðtÞ, IðtÞ, 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/6 y , γ = 1/3d = ð365/3Þ y −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.

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.
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, 2 Computational and Mathematical Methods 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).
Then, the sequential period doublings appear in the bifurcation diagram, Figure 2, 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  3 Computational and Mathematical Methods 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.

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 Figure 3(a). After a very large outbreak, we observe rather small outbreaks next, a large I n followed by a small I n+1 .
For intermediate small I n , we observe rather large I n+1 . But for very small I n , also, the I n+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.
Though such discrete maps generated from 3dimensional ODE flows have to be two dimensional, as seen in Figure 3(a), the described sequence of I n+1 following I n 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 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]. The period doubling bifurcations cumulate into a fuzzy clowd of points in the bifurcation diagram, Figure 4(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" 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 around a ≈ 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.

Computational and Mathematical Methods
Namely, we determine subsequent bifurcation points from zooming into the bifurcation diagram in smaller and smaller regions of the parameter a, see, e.g., Figure 5 for the final step, and observe a c1 = e 2 = 7:389056099, which can be calculated analytically, then a c2 = 12:509, giving a difference of Δa 2 ≔ a c2 − a c1 = 5:12.
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.

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

Computational and Mathematical Methods
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, 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 fx n g 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).
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., x n+2 = f ð f ðx n ÞÞ ≕ f ð2Þ ðx n Þ, 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 x n+4 = f ð4Þ ðx n Þ, 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 x n = 0:5, compare Figures 6(f) and 6(b). Now, the original map Figure 6(f) can be also compared with the second iterate Figure 6(d), when we flip the inner part around x n = 0:5 and zoom in or mathematically rescale x-axis and y-axis. This rescaling can be given mathematically in the following way: 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 − f ð2Þ 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, x n+1 = 1 − μx 2 n ≕ f μ ðx n Þ is generally used with the approximate self-similarity with scale factor again called α. This defines an operator T which maps 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 gðxÞ for period doubling going to infinity fulfilling which gives explicitly the so-called Cvitanović equation for the universal constant α = 2:5029 and universal function gðxÞ. 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 T ½gðxÞ = gðxÞ, 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 of T around the renormalization fixed point T g = g by for small perturbation hðxÞ, resulting explicitly in and obtain the eigenvalue/eigenfunction problem in function space with identified at the onset of chaos, after infinitely many period doubling bifurcations, at parameter value a with value a ∞ = 3:5699456.
The fixed point equation T ½gðxÞ = gðxÞ then can be solved via popynomial ansatz of gðxÞ giving also the value for α and the linearization L½hðxÞ = δ F · hðxÞ 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].

Practical Solution of the Functional Equations.
We start with the polynomial Ansatz for the universal function g 6 Computational and Mathematical Methods quadratically via gðxÞ = g 0 + g 1 x + g 2 x 2 + Oðx 3 Þ to be inserted into T½gðxÞ = gðxÞ = −αgðgð−x/αÞÞ. This gives via coefficient comparison in orders of x the results g 1 = 0, g 0 = 1 and g 2 = −α/2, and from this compared to the higher order polynomial results of 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, hðxÞ is inserted into L½hðxÞ = δ F · hðxÞ which already gives in zeroth order compared to higher order result of 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.

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.

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 The master equation can be written in a generic form using densities x 1 ≔ S/N and x 2 ≔ I/N instead of the total numbers of individuals in the population classes S and I, hence state vector x ≔ ðx 1 , x 2 Þ tr , as with n = 3 different transitions and small deviation from state x as Δx j ≔ ð1/NÞ · r j . For the master equation in densities x 1 ≔ S/N and x 2 ≔ I/N, we have now transitions w j ðxÞ and its vectors r j given by 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].

Fokker-Planck Approximation of the SIR Model.
From the master equation in densities, Equation (24) for the SIR system, we use Taylor's expansion giving to second order in 1/N a Fokker-Planck equation with or in a different notation 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].

Conclusions
We have investigated models able to describe respiratory disease outbreak dynamics, in cases relevant for influenzalike 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 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. models in completely different parameter regions, describing childhood diseases [31], but the interplay between such coexisting 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.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.