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.
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: and its corresponding dynamics of probabilities is discussed below. The population of size is divided into susceptible individuals , infected , and recovered , 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 and finally,
This stationary state is a fixed point in the state space of the variables and .Adding the seasonal forcing 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 , , and . 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 , , , and and varying . For an analysis based on the time available influenza data from France, we took as analyzed in detail in an earlier publication [2]. Here, we consider from the nonforced system, i.e., , and no import a gradual increase of the seasonality, and observe already for very small seasonality bifurcations and complex dynamics, much below .
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., , 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., , 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.
(a)
(b)
(c)
Then, already in a very close by seasonality of , 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 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 , 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 , , and . From these first values, we can calculate a “speed of bifurcation accumulation” . 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 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 Figure 3(a). After a very large outbreak, we observe rather small outbreaks next, a large followed by a small . For intermediate small , we observe rather large . But for very small , also, the 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.
(a)
(b)
Though such discrete maps generated from 3-dimensional ODE flows have to be two dimensional, as seen in Figure 3(a), the described sequence of following can be roughly represented by a simple functional form of . 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 , see Figure 4(a), giving the bifurcation diagram of the Ricker map for, e.g., [4, 9].
(a)
(b)
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 to give a Lyapunov spectrum, see Figure 4(b). We use the same range of the parameter 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 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 parameter , see, e.g., Figure 5 for the final step, and observe , which can be calculated analytically, then , giving a difference of .
Continuing this analysis further, we obtain , , then , , and , , and finally , . Hence, for the speed of bifurcations, we obtain , , , and . Hence, we move towards a , which gives for subsequent bifurcation points , ones , measuring the convergence speed into chaos at .
Remember that we now found values of the speed of bifurcations close to 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 . 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 towards .
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 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, 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 , the universal Feigenbaum constant , 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., , and observe that the trajectory converges to a fixed point, but for a value of, e.g., keeps oscillating in a period 2, see Figure 6(e).
(a)
(b)
(c)
(d)
(e)
(f)
Then, for, e.g., , we find a period 4 cycle in the Feigenbaum map, Figure 6(c). If we iterate the map, i.e., , the period 4 becomes a period 2 limit cycle again, see Figure 6(d). And for, e.g., , 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 , 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 , 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 and zoom in or mathematically rescale -axis and -axis. This rescaling can be given mathematically in the following way: with a scaling factor (or here simply since it cannot be confused with the parameter in the SIR models) and flipping via in the -axis and in the -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, is generally used with the approximate self-similarity with scale factor again called . This defines an operator which maps into 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, has a fixed point function for period doubling going to infinity fulfilling which gives explicitly the so-called Cvitanović equation for the universal constant and universal function . The constant is the first Feigenbaum constant, while however we are here more interested in the second Feigenbaum constant . The renormalization theory now gives an attractive high-dimensional manifold, attracting towards the fixed point , and one unstable direction in function space with dominating eigenvalue, which turns out to be we are looking for [16].
Hence, we investigate the linearization of around the renormalization fixed point by for small perturbation , 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 with value .
The fixed point equation then can be solved via popynomial ansatz of giving also the value for and the linearization giving in addition , 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 , 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 function quadratically via to be inserted into . This gives via coefficient comparison in orders of the results , and , 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, is inserted into 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 , in terms of bifurcation accumulation 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 . 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
The master equation can be written in a generic form using densities and instead of the total numbers of individuals in the population classes and , hence state vector , as with different transitions and small deviation from state as .
For the master equation in densities and , we have now transitions and its vectors given by
We can perform simulations of the master equation directly using the Gillespie algorithm [18, 19], observing that for large population sizes 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 giving to second order in a Fokker-Planck equation with or in a different notation using simply a quadratic form here with
The Fokker-Planck equation gives a stochastic differential equation system with , and in the SIR model, the two-dimensional Gaussian normal noise vector as and using matrix square root from eigenvalue-eigenvector decomposition as 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 being a quadratic matrix and allow a matrix with , such that the expressions in can become quite simple and fast to compute [22]. Explicitly, we have for the expression in the SIR model given as , giving again as required.
In Figure 7, we show simulations for the seasonal SIR model with the previously stated parameters and high seasonality , for the mean field ODE system, the master equation simulated via the Gillespie algorithm, and the stochastic differential equation systems using the quadratic matrix and the nonquadratic 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.
(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 the number of reactions by Poisson distributions, assuming little changes in the transitions during the interval , the so called -leap algorithms, which tend for increased from the Poisson distributions again to the Gaussian distributions of the stochastic differential equation systems, as given in Equation (32) with the -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 method to evaluate waiting times after the event of some characteristic reactions, eventually strengthening the condition [26]. Here, the -distribution as continuous time limit of discrete time Erlang -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 , 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 during the simulations. Basically, in the simplest case of a one-dimensional Langevin equation of the typ , the state-dependent noise term is evaluated in the interval as with , giving in its simplest version for the Ito version and for the Startonovich version [28], and adding to the drift the correction to Ito type, namely following obviously from the Fokker-Planck formulation [9, 20], gives the Milshtein method, either kept implicit to desired order in 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We thank Stephan Thomae and Heiner Müller-Krumbhaar for introducing some of us to universality concepts and especially to the Feigenbaum renormalization and also thank Keith Briggs for further discussions on numerical aspects of the Feigenbaum constants and Pedrac Cvitanović for discussions on renormalization, both in stochastic processes and chaotic dynamics. We would like to thank the Basque Health Department and the Basque Health Service (Osakidetza) for the first data on influenza notifications and for firsthand information and rapid discussions of the ongoing efforts to combat the COVID-19 crisis. We thank the huge efforts of the whole COVID-19 Basque Modelling Task Force (BMTF), especially Adolfo Morais Ezquerro, Vice Minister of Universities and Research of the Basque Government for starting the BMTF initiative and for fruitful discussions. M. A. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 792494. This research is also supported by the Basque Government through the BERC 2018-2021 program and by the Spanish Ministry of Sciences, Innovation and Universities: BCAM Severo Ochoa accreditation SEV-2017-0718. N.St. thanks the Mathematics Department of Trento University for the opportunity and support as a guest lecturer in biomathematics, during which part of the present work has been conceived and started to be developed.