How can contemporary climate research help to understand epidemic dynamics? -- Ensemble approach and snapshot attractors
ssnapshot attractors in SIR
How can contemporary climate research help to understand epidemicdynamics? – Ensemble approach and snapshot attractors
T. Kovács a) Institute of Theoretical Physics, Eötvös University, Pázmány P. s. 1A, H-1117 Budapest,Hungary. (Dated: 13 August 2020)
Standard epidemic models based on compartmental differential equations are investigated under continuous parameterchange as external forcing. We show that seasonal modulation of the contact parameter superimposed a monotonicdecay needs a different description than that of the standard chaotic dynamics. The concept of snapshot attractors andtheir natural probability distribution has been adopted from the field of the latest climate-change-research to show theimportance of transient effect and ensemble interpretation of disease spread. After presenting the extended bifurcationdiagram of measles, the temporal change of the phase space structure is investigated. By defining statistical measuresover the ensemble, we can interpret the internal variability of the epidemic as the onset of complex dynamics evenfor those values of contact parameter where regular behavior is expected. We argue that anomalous outbreaks ofinfectious class cannot die out until transient chaos is presented for various parameters. More important, that thisfact becomes visible by using of ensemble approach rather than single trajectory representation. These findings areapplicable generally in nonlinear dynamical systems such as standard epidemic models regardless of parameter values.PACS numbers: SIR model, seasonal forcing, ensembles, snapshot attractor, transient chaos
I. BACKGROUND
The effort to estimate short and long-term behavior, systemparameters, possible control strategies, and also social impactof COVID-19 pandemic stimulated an explosion of diseasemodelling and scientifically demanding examinations frommany different point of view made by various communitymembers and research groups around the globe. Althoughthe present situation is rather new, we have many examplesfrom the history that helped to lay down basic mathematicalrules of disease propagation. The celebrated SIR (suscepti-ble - infectious - recovered) model splits the population intothree disjoint parts and deals with the number of individualsin these sub-populations as time goes on.The SIR model and its variants [SI, SIS, SEIR, SEIRS,RAS] show qualitatively similar dynamics and are in goodagreement with observations. In a homogeneous environmentthese models possess a globally stable fixed point attractoras a disease equilibrium . Omitting the stochastic nature ofreal world disease spread , the deterministic SIR-like mod-els present diverse and reach dynamics. The non-linearity ofthe model and also the time dependent internal forcing canbe considered as a source of complexity.Sufficiently large seasonal forcing or different mixing ratesbetween sub-populations can cause large oscillations and alsoperiod doubling cascades . It has also been demon-strated that for certain values of system parameters the non-autonomous models of recurrent epidemics (measles, mumps,rubella, H5N1 avian influenza) show chaotic behavior .Traditional SIR-like epidemic models are dissipative nonlin-ear low-dimensional systems with constant, periodic, quasi-periodic , or term-time external forcing whose dynamics a) Electronic mail: [email protected] is governed by (chaotic) attractors in phase space. Further-more, even if the long-term dynamics is regular and the finalstate of the system is a fixed-point attractor the route to thiscondition might be rather complex. Many studies examinedthe role of finite-time irregularity in ecological models as well as in epidemic dynamics concluding the rele-vance of transient behavior.It is known from dynamical systems theory that chaotic at-tractors always can be characterized by the natural measurecorresponding to the distribution of possible states in phasespace . It turned out, however, that traditional numericalmethods, such as monitoring a single chaotic long-term tra-jectory, fail in case of arbitrary forcing (like smoothly vary-ing parameters). Obviously, if one wants to model such asystem (e.g. the changing climate with shifting atmosphericCO concentration), this shortcoming has to be overcome.The mathematical concept of snapshot attractors , knownfor many years in (theoretical and experimental) dynamicalsystems community , fulfills entirely our wish.Briefly, snapshot attractors are time-dependent objects inphase space of non-autonomous dynamical systems. Thus theshape of snapshot attractors is changing in time while theirfractal dimension might even remain constant . Furthermore,obtaining one single trajectory in a system with arbitrary driv-ing force, does not provide the same result as the ensembleapproach (many trajectories emerging from slightly differentinitial conditions) in the same system at a given time instant.This effect is the consequence of the fact that ergodicity is notsatisfied as the system is driven aperiodically . This con-clusion generated the recent opinion that the changing cli-mate should be scrutinized by ensemble approach ( parallelclimate realizations ) rather than averaging a single long-termtime series .In this work a seasonally forced deterministic epidemicmodel with monotonically changing contact rate (due to, forinstance, vaccination or restricting the movement of the popu- a r X i v : . [ q - b i o . P E ] A ug napshot attractors in SIR 2lation) is presented. We propose that the statement by Ref. 43 ”climate change can be seen as the evolution of snapshot at-tractors” also holds for disease spread dynamics with contin-uously changing contact parameter.In Section II the epidemic model is defined. After that,the results about stationary (Sec. III) and changing epidemic(Sec. IV) are presented. Section V is devoted to conclusions. II. STANDARD EPIDEMIC MODEL
Compartmental disease models describe the number of in-dividuals in a population regarding their disease status: sus-ceptible ( S ), infectious ( I ), or recovered ( R ). Although, thesemodels involve many simplifications (such as the progressionof infection or difference in response of individuals) they per-formed well in real world epidemic situations . There aretwo major groups of epidemic models: the SIR-like clusterthat characterizes lifelong immunity (e.g. measles, whoop-ing cough), and the SIS-like class (containing mostly sexuallytransmitted diseases) which portrays repeated infections .Here we study the SEIR equations that involves a fourthgroup in addition to previous ones. More concretely, we as-sume that an individual enters the population at birth as sus-ceptible and leaves it by death. A susceptible becomes ex-posed ( E ) when contacting one or more persons, called infec-tive(s), that can transmit the disease. In an incubation periodthe exposed individuals are infectious but are not yet infec-tious. After this term they become infective and become laterimmune or recovered.The mathematical model associated with above descriptionreads as follows d S d t = m − bSI − mS , d E d t = bSI − ( m + a ) E , d I d t = aE − ( m + g ) I , d R d t = gI − mR (1)with the following notations and assumptions : (i) S , E , I and R are smooth functions of time and the size of the whole pop-ulation remains unchanged, S + E + I + R =
1; (ii) there areequal birth and death rates ( m ); (iii) the probability that anexposed individual remains in this class for a time period τ after the first contact is exp ( − a τ ) , where 1 / a is the mean (orcharacteristic) latent period; (iv) similarly 1 / g gives the meaninfectious period, the time period that an individual spends asinfectious before recovery; (v) immunity is permanent and re-covered individuals do not re-enter to susceptible class. Thecontact rate b ( t ) , the average number of susceptibles con-tacted per a single infective per unit time, is the origin of thespread of the disease. In the case of annual periodicity b ( t ) = b ( t )[ + b cos ( π t )] . (2)and t is measured in units of years. When the system parameters are kept constant( m , a , g , b , b = b > g ). In this case the dynamics can be further characterizedby the expected number of secondary cases caused by aninfectious individual in the susceptible population. This isknown as the basic reproductive ratio, and can be expresseddenoted here by R = ba / [( m + a )( m + g )] . R essen-tially determines whether a disease can ( R > , endemicequilibrium) or cannot ( R < , infection-free) persist in apopulation. Provided that latent and infectious periods areshort, i.e. a , g (cid:29) m , one can use the approximation R = b / g . There are observations, for example in childhood diseasesor avian influenza, that do not show damped characteristicsrather (ir)regular cycles instead. This phenomenon can belinked to seasonal variation in contact rate b ( t ) or in recruit-ment rate g . It has already been noticed that in case of peri-odic forcing, b ( t ) = b = const. in Eq. (2), chaos emergesin epidemic time series . Similarly to the Lorentz84model wherein the solar irradiation induces seasonal ef-fect, the SEIR model with variable contact rate is a non-autonomous low-order system with well-defined chaotic at-tractor in the phase space at certain parameter values. Thelow dimension of the phase space allows to visualize its pat-tern comfortably.As a result of the fixed population size, point (i) above, oneof the four equations (traditionally the equation for R ) canbe omitted. In addition, the infectious and exposed groupsturn out to be linearly related to the first order at least,see Figure S1 in supplementary material. Consequently, the S − I phase portrait represents the state space texture accu-rately. Plotting one trajectory in the S − I plane of the phasespace one would observe a coil design due to the explicit timedependence of the model, b (cid:54) =
0. A conventional methodto make a periodically forced system autonomous is to take”pictures” about the phase space with the same frequencyas the excitation acts, i.e. making a stroboscopic map .By choosing an initial time t and purely sinusoidal force inEq. 2 ( b ( t ) = const.) one can interpret the filamentary shapeof the chaotic attractor by defining the stroboscopic map t mod T = . Therefore, the periodic variation in Eq. (2) yieldsam annually stationary epidemic with steady attractor.Eq. (2) assumes a constant average contact rate b . We havean annual period starting with its maximum at t = . Inphysical or biological context this means the largest contactrate possibly due to school term, vacation and holidays, sea-sonal breeding pattern in a seabird colony, etc. Eq. (1) withthe following parameters, corresponding to measles , showsirregular dynamics: m = .
02 year − , a = .
84 year − , g =
100 year − , b = − , b = .
28 year − . The aboveparameters lead to R ≈ . Numerical calculations show thatfor lower mean values of the contact rate periodic fixed-pointattractors exist. For instance, biannual cycles at b ≈ . For more detailed asymptotic dynamics see the bifurcation di-agram in Figure 1 (blue curve).Until now b was kept as a constant. Following the climatechange methodology monotonic variation in b ( t ) through b ( t ) is established. Equation (3) specifies the temporal de-napshot attractors in SIR 3cay of the mean contact rate b ( t ) = (cid:40) t ≤ t st e − α ( t − t st ) +
90 if t > t st . (3)Here t st represents the time until the epidemic is stationary,i.e. b = const., this is set to be 250 years after initialization.This amount of time seems to be enough since the conver-gence of trajectories to the attractors is much shorter, t c ≈ α is the decay rate of b ( t ) . In this study α takes three values 0,04 yr − , − , and 0.004 yr − , and b always starts at 1800, i.e. from the chaotic regime. Eq. (3)implies that lim t → ∞ b ( t ) = . This means that the criticalvalue of mean contact rate b =
100 (or in terms of basic re-production rate R ≈
1) is reached at different times after t st ,according to α . The main motivation of this study is based on Eq. (3). Con-sidering, for example, the time series of childhood measlesin London after World-War II (between 1945 and 1990), thevaccination program (starting in 1968) changed dramaticallythe number of cases . Beside the medical treatment, ifavailable, other artificially forced decline in contact rate (suchas gradual lockdown prescribed by administration) can be im-itated by Eq. (3) resulting in a changing epidemic. This secu-lar variation of parameter b gives a clear analogy to climatechange models and the framework applied in this field. III. STATIONARY EPIDEMICA. Parameter dependence
The bifurcation diagram basically reveals the long-term dy-namics, excluded the initial transients, in a given parameterrange. Classically, one single trajectory sampled by strobo-scopic map serves the blue shape of bifurcation diagram inFigure 1. In case of the stationary or stationary epidemic( b = const.) complex dynamics arises for certain parame-ter values in the SEIR model, b (cid:38) .The same applies for the phase space pattern, too. Fig. 2 de-picts the chaotic attractor obtained after T = S − I phase por-trait is pictured in Figure 2. The chaotic attractor ( b = T = S − I plane. FIG. 1. Bifurcation diagram of the SEIR model with stationary b .The variable S is plotted on the same year of the day starting at t = . The classical diagram (blue curve) is obtained by integrating 2000initial conditions for t = b ≈ . Theend-points of the same trajectories are plotted after 500yr iterationindicating long-lived transients (green). Clearly finite time chaoticmotion unfolds the extra structure of bifurcation diagram. Upperx-axis evaluates the basic reproductive ratio derived from system pa-rameters. R = R < b < R of common diseases forcomparison only. (Ebola ∗ : based on 2014 Ebola outbreak.) Transient chaos, complex behavior on finite-time scales ,manifests before a trajectory comes to an attractor, both indissipative as well as in conservative systems. In general, theattractor can be a simple object in the phase space, for in-stance, a fixed-point or a limit cycle. Following the time evo-lution of the ensemble before it approaches the attractor, onecan capture other ingredients of the bifurcation diagram byconsidering the initial transients. To do this, the uniformlydistributed 2000 trajectories are integrated for 500 years andthe corresponding state space positions are stored. This fea-ture becomes visible in the green vague domain that extendsalong the blue curve. Clearly, transient chaos has a significantcontribution to the dynamics for b (cid:38) . B. Ensemble view
The initial transients are usually thrown away in long-termdynamical analysis. However, the complexity might appearjust in this phase of motion as indicated by the bifurcationdiagram in Figure 1, already in the stationary case. Illus-trating the role of finite time chaotic behavior in the SEIRnapshot attractors in SIR 4
FIG. 2. The stationary chaotic attractor in the S − I phase plane.The ensemble is plotted at three different time instants: t = , , (cid:104) S (cid:105) and (cid:104) I (cid:105) values from Fig. 3 (middle panel) while theellipse indicates the associated standard deviations (bottom panel),respectively. Note that the green circle is a weighted average (a kindof ’barycentre’) of the trajectories along the attractor rather than itsgeometrical focus. epidemic model we track the evolution of 1600 different ini-tial conditions distributed uniformly in a cube S = [ . ] , E = [
0; 0 . ] , I = [
0; 0 . ] , .Fig. 3 demonstrates how the individual members of theensemble reach the attractor at different time instants. Forthose values of contact parameter when the epidemic showsbiannual ( b = b = (cid:104) A ( t ) (cid:105) = / N ∑ Ni = A i ( t ) (where A i ( t ) denotes an observablecorresponding to the i th member at time t) designates how fastthe initial irregularity terminates and also adverts the averagetransient time. One can observe the similar characteristics forboth the (cid:104) S (cid:105) and (cid:104) I (cid:105) curves, fixed-point and chaotic attractors(top and middle panels), respectively.In the stationary epidemic model one might expect that aftersome time the extent of the attractor remains constant suggest-ing that all of the individual trajectories in the ensemble havereached it. The standard deviation, σ A = ( (cid:104) A (cid:105) − (cid:104) A (cid:105) ) / , (4)refers to the size of the attractor in the direction of A . Figure 3(bottom panel) depicts σ S and σ I in case of b = . Theinitial small values of standard deviations show that the en-semble moves together at the beginning of integration. Then,the members spread out in the phase space ( ≈ −
50 years)and after a while start to approach the attractor. This time islonger for parameters b = σ S is smaller. FIG. 3. Ensemble view of stationary dynamics. Top and middle:the mean of S (blue) and I (red) variables vs. time. Some of theindividual trajectories are also shown sampled by stroboscopic map(only the S co-ordinate). The individual trajectories arrives sooneror later the biannual fixed-point attractor (see for example the cyanparallel dots for t >
230 in top panel). Bottom: standard deviationsof the same variables as in the middle for b = Stationary epidemic with constant σ A essentially means thatafter the transients the phase space structure does not changein time. This can be visualized by using of ensemble ap-proach. Without loss of generality we always start our sim-ulation at t = b ( t ) ) and take the next pictureabout the ensemble at t mod T = natural distribution of a chaotic attractor comes fromthe fact that the dots visit certain parts of the filamentary struc-ture more frequently than others. Moreover, this distributionis stationary and does not depend on initial conditions. Insteadof investigating a 3D frequency diagram (see supplementaryFigure S2), we plot the projection I of trajectories wander-ing on the chaotic attractor in Figure 2. To obtain Figure 4,a fine grid is defined in the S − I phase plane and the num-ber of points in each cell is plotted against the variable I . Thehistogram consists of four different ensembles. Ensemble 1,2,napshot attractors in SIR 5
FIG. 4. Natural distribution on the chaotic attractor ( b = I variable after t = and 3 covers the same volume in the phase space but involvesslightly different initial state vectors. Ensemble 4 lays out ini-tial conditions from other part of phase space. The distributionclearly shows the same pattern for all four ensembles. IV. EPIDEMIC CHANGE
In a changing epidemic the parameter(s) of the systemis(are) varying as time goes on. One might think that, say,a decreasing of the contact parameter b means to walk alongthe bifurcation diagram slowly from right to left, and terminat-ing at a safe destination with R < . In this section we willshow that this idea is fairly naive due to the internal variabilityand the transient effects in the epidemic. To capture the dy-namics properly in this scenario, the ensemble approach andthe concept of snapshot attractor is desirable.
A. Snapshot attractor geometry and natural distribution
According to the common measles parameters we start theswitch-off process, called also epidemic change, from thechaotic attractor ( b = reveal that starting from the chaotic attrac-tor well after that the trajectories reached it, the evolution ofdifferent states might be rather diverse. That is, transient dy-namics becomes important again while the parameter changeappears in the system. In Figure 5, four trajectories are se-lected from the chaotic attractor and their time dependenceis followed. It can be immediately seen that different colorsreach the fixed-point attractor at different times. For instance,in case of α = .
004 (bottom curve), first the blue, then theyellow, red, and green curves arrive at the fixed-point attrac-tor. One can pick up other trajectories that will have longeror shorter oscillations. Furthermore, the previous order ofcolors, i.e. transient times, might change with α , as demon-strated by the middle curve. Therefore, we can point out, thatit is worth investigating several trajectories simultaneously in-stead of following individuals. However, after the parametershift is switched on, the shape of the chaotic attractor starts to FIG. 5. Four trajectories plotted in blue, green, red, and yellow wan-dering on the chaotic attractor and their fate (in S variable) after pa-rameter change sets in at t =
250 yrs (vertical line). The two scenar-ios α = . , .
04 are shifted by 0.2 units, respectively, for bettervisualization. The very peculiar time evolution of individual trajec-tories for different α s suggests to use ensembles. vary and forms a time-dependent set called snapshot attrac-tor. Thus, snapshot attractor is an object in the phase planethat contains the whole ensemble at a given time instant. Wenote that this mutation is not a result of the particular choiceof initial phase (day of the year) of the mapping rather thanthe change of parameter. The geometrical alteration of thesnapshot attractor is followed by the change of distribution onit. Nevertheless, its form and the distribution is independentof the choice of the ensemble and the onset ( t st ) of parameterdecay.Figure 6 exhibits snapshot attractors drawn by N = ini-tial conditions during the epidemic change. The mean contactrate decreases from top to bottom and left to right as speci-fied on each panel. Thus, every single plot corresponds to adifferent time instant of simulation. The top left panel, takenat T = t < t st = b is well below of its original value (1800) as well as of 1770where large extended chaotic attractors are formed in the bi-furcation diagram (Fig. 1 inset). In other words, for those con-tact rates ( b = , α = .
01 one can obtainsimilar alteration for different switch-off rates too. Slowerscenario (e.g. α = .
004 in Eq. (3)) allows the attractorto keep its filamentary shape and maintain chaotic dynamicslonger. From other perspective, the same parameter value b is reached later while α is smaller. The opposite is true for afaster parameter change, say α = . . It can also be shownthat approximately the same pattern belongs to the same con-tact value regardless the rate of change. That is, the fasterthe contact rate decay, the less pronounced the complexity inphase space.napshot attractors in SIR 6
FIG. 6. The evolution of the snapshot attractor with α = . . From left to right and top to bottom the ensemble is taken at t = , , ,
618 yrs. The first picture still corresponds to stationary dynamics (before t st =
250 yr) with the original b = . Filamentary pattern persist for lower b values, too, where no stationary chaotic attractor exists in phase space according to the bifurcationdiagram (Fig. 1). Still to be noted the physical extension of the snapshot attractor that is also time-dependent and their size varies significantlyin both directions. Compare the scale of axes in different panels. The green circles and the red ellipses denote the same quantities as inFigure 2. Slim panels to the right of the S − I phase planes display the natural distribution projected onto the I variables. For b =
133 mostof the phase points accumulate around the fixed-point attractor ( I ≈ − ), therefore, no histogram is presented. Initial conditions cover thefollowing state space volume: S = [ . ] , E = [ . ] , I = [ . ] . Other interesting feature is that although the contact rate b decreases monotonically, the size of the snapshot attractorsmight increase, Figure 6 top right panel. This fact yields thatthe domain accessible by the dynamics can be larger. Thusthe possible ( S , I ) pairs extend to larger domain in phase spaceeven for smaller contact parameter. Furthermore, the average(green circles) and the standard deviation (red ellipses) alsochange in time. And this temporal behaviour cannot obtainedfrom the classical view, only by using snapshot attractors.To understand this phenomenon we should recall the con-cept of chaotic saddle. This non-attracting set with its sta-ble and unstable manifolds in phase space is responsible forthe finite time chaotic behavior . To construct the saddlenumerically, we define two holes in the phase space (blackrectangles in Figure 7). Then, a large number of initial con-ditions ( N = · ) are distributed uniformly in the region S = [ .
04; 0 . ] , E = [
0; 0 . ] , I = [ . . ] and thetrajectories a integrated for T =
30 years. Only those trajecto-ries are kept that do not enter the holes during the integrationtime. The initial conditions belonging to these trajectoriesdraw the stable manifold of the saddle (red dots in Fig. 7),while those that are just before leave plot the unstable man-ifold (green dots.) The saddle itself is the intersection of its manifolds, not shown here. The average lifetime of chaos (theinverse of escape rate, κ ) can be estimated by calculating thetime distribution of the non-escaped trajectories.It is also known from the theory of transient chaos that thesaddle’s manifolds have filamentary design just like the snap-shot attractors in Fig. 6. Analogously to chaotic attractor, sta-tionary dynamics with constant driving amplitude defines astationary chaotic saddle related to certain parameter values ofthe system. Due to the continuous adjustment of the contactrate, b ( t ) , in the changing epidemic model, the trajectroriesdo not have time to reach the attractor belonging to a given b value. Consequently, a time dependent chaotic saddle is con-sidered, whose unstable manifold approximates the snapshotattractor . This stationary saddle persist for very low b val-ues and its unstable manifold controls the finite time complexepidemic dynamics, Figure 7.We emphasize at this point that the natural distribution as-sociated to the snapshot attractor also changes in time follow-ing the geometrical reorganization of the phase space pattern(histogram visualization in Fig. 6).napshot attractors in SIR 7 FIG. 7. The snapshot attractor (blue) and stable/unstable (red/green)manifolds of the stationary chaotic saddle at b = . Discon-tinuities along the unstable manifold reflect to the escape condi-tions of the numerical scheme. The long-term dynamics of the sta-tionary model is visualized by the four tiny yellow circles around ( S , I ) ≈ ( . , . ) as fixed-point attractors. The snapshot at-tractor coincides with the top right phase portrait in Figure 6 B. Parallel epidemic realizations
Similarly to climate research we can define the concept of parallel epidemic realizations.
In standard disease modelslike Eq. (1) this can be done by using a large number of initialconditions as an ensemble in phase space and following theirevolution as discussed in Section III A. This picture can beimagined as many copies of the epidemics obeying the samephysical laws and being affected by the same time-dependentforcing .As we have seen before, in case of a changing epidemic,the time-dependent forcing has an impact on the natural mea-sure of the snapshot attractors. This fact results in a temporalchange of the average values as well as internal variability.To quantify the internal variability of the system (1) statisti-cal measures over the ensemble should be proposed like inthe case of stationary epidemic. The variance or the standarddeviation σ A in Eq. (4) of the ensemble characterize the fluc-tuations around the averages indicating the inherent internalvariability.Top panels of Figure 8 show the ensemble standard devia-tion σ S (blue) and σ I for different decay rates α . The paral-lel epidemic realizations contain 1600 numerically integratedtrajectories sampled at integer multiplies of years. The con-vergence time t c , the time needed of the ensemble to reachand spread along the attractor, is found to be t c = σ S and σ I are constant before t st = b = − . Right afterthis the contact rate starts to decrease ( t > σ S ( t ) increases, after the maximum the trend follows nearly b ( t ) , it decaysto zero illustrating that the size of the attractor shrinks andasymptotically reaches the neighbourhood of a regular fixed-point attractor, ( S , I ) = ( , ) , as expected for b ≈
90 ( R = . α , the more regular, i.e less filamentary,the phase space geometry at the same time, see bottom panelsin Fig. 8.The shape of the standard deviation can be explained bytransient chaos that occurs during the parameter change. Thesize of the green shaded band around the main feature in thebifurcation diagram (Fig. 1) already indicates that the size ofthe phase space region filled with transient chaos increases as b reduces. To demonstrate this we call the attention to thehorizontal dimension of the snapshot attractors in Figure 6.The smoothly decreasing profile of the standard deviation af-ter reaching the maximum refers to the existence of transienteffect up to very small parameter values.Different maximum values of the σ S ( t ) curves indicate anon-trivial relation between the internal variability and chang-ing rate of b . The largest value corresponds to α = .
01 (topmiddle panel), the other two scenarios ( α = .
04 and 0.004)show roughly the same maximum size of the snapshot attrac-tor, albeit, at different time instants. Bottom panels present themaximal physical extension of the snapshot attractors versustime. These plots also supports our observations that the sizeof the attractor first increases and the shrinks to be a fixed-point attractor.One possible explanation of this property is the relative”coupling” between the timescales, that is, decay of forcingdefined by the exponent α and the average lifetime of chaos instationary dynamics at specific parameter b . From a geomet-rical point of view, how close the snapshot attractor evolvesto the unstable manifold of the stationary chaotic saddle atparticular contact parameter. A detailed exploration of thisfeature is, however, postponed to a future study.
V. FINAL REMARKS
Parallel climate realizations, an effective and new frame-work in climate research, has been adopted to an epidemicmodel to explore the fading of the complexity due to the sys-tematic switch off the driving mechanism. The mathematicalconcept of snapshot attractors and their natural distributiondemonstrates that single time series analysis is not capableto reflect the complex dynamics of a changing epidemic. In-stead of monitoring isolated events, the ensemble view of tra-jectories – parallel epidemic realizations – and its statisticaldescription is desirable.The temporal change of the attractor geometry and the dis-tribution on it reveals the internal variability of the dynamics.The extension of the bifurcation diagram indicates the im-portance of transient chaos in the stationary dynamics as wellas during the continuous parameter shift. No matter whetherwe start from a stationary state of the long-term dynamics, theswitch-off process activates the hidden parts of the bifurca-tion diagram. Thus, the underlying non-attracting object, thechaotic saddle, or more precisely its unstable manifold, orga-napshot attractors in SIR 8
FIG. 8. Statistical measures of SEIR-ensembles. Top: The standard deviation of variable S and I for various decay processes. The epidemicchange starts at t st =
250 yrs well after the convergence time ( t c ≈
50 yrs). Vertical dashed lines mark the time instants corresponding to phasespace portraits in Figure 6. Bottom: The linear size of snapshot attractors in both variables indicate the role of transient chaos and internalvariability of dynamics depending on the decay rate α . nizes the system’s evolution.The dissipative relaxation timescale, the inverse of thephase space contraction rate based on the divergence of thevector field of system (1), is extremely low, ∼ − α − ≈ −
250 yrs, κ − ≈ −
35 yrs, re-spectively. This latter is obtained a few particular contact rate.This implies that the switching off process, with parameter α used in this study, is not quasistatic . In other words, there isnot enough time for trajectories to reach the stationary attrac-tor, either it is chaotic or regular, due to the nonstop parametervariation. Therefore, the dynamics of the switching off pro-cess shows much more complexity than that of the stationaryscenario.In the spirit of the investigated SIR-like classical epidemicmodel (1) with parameter shift of the contact rate b ( t ) weemphasize the importance of the finite time irregular dynam-ics. Our results reinforce the ”hidden” complex transients,and also the importance of time-dependent snapshot attractorsnot just for an epidemic with large reproductive ratio, suchas measles and chickenpox, but for those with lower values( R (cid:38)
1) too.The main conclusion of this study sheds light on the im-portance of the ensemble view and parallel realizations in epi-demic dynamics. Fig. 6 illustrates the general features oneexpects in the (long-term) prediction of any epidemics. Pre-dictions based on individual simulations are not reliable since,due to the chaotic nature of the dynamics, they can lead tomany possible results, to any point of the snapshot attractorbelonging to the time instant of the prediction. The ensembleapproach is able to treat all possible outcomes as a whole, andpredicts even the probability of the different permitted epi- demic states. When convenient, statistical moments of thisdistribution can be determined, e.g. the average (providingthe most typical epidemic outcome), or the variance (charac-terizing how broad the distribution of the permitted states is),but higher order moments might also be useful. Anyhow, onecan thus follow how the statistical prediction changes in time.
ACKNOWLEDGMENTS
We benefited from useful discussions with T. Tél. Thiswork was supported by the NKFIH Hungarian GrantsK125171. The support of Bolyai Research Fellowship andÚNKP-19-4 New National Excellence Program of Ministryfor Innovation and Technology is also acknowledged. See the pre-filtered search on http://arxiv.org andhttps://connect.biorxiv.org/relate/content/181. W. . Kermack and A. G. McKendrick, “A contribution to the mathematicaltheory of epidemics,” Proceedings of the Royal Society A , 700–721(1927). M. J. Keeling, P. Rohani, and B. T. Grenfell, “Seasonally forced diseasedynamics explored as switching between attractors,” Physica D NonlinearPhenomena , 317–335 (2001). H. C. Tuckwell and R. J. Williams, “Some properties of a simple stochasticepidemic model of sir type,” Mathematical Biosciences , 76–97 (2007). L. J. S. Allen, “An introduction to stochastic epidemic models,” in
Mathe-matical Epidemiology , edited by F. Brauer, P. van den Driessche, and J. Wu(Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 81–130. L. F. Olsen and W. M. Schaffer, “Chaos Versus Noisy Periodicity: Alterna-tive Hypotheses for Childhood Epidemics,” Science , 499–504 (1990). P. J. Witbooi, “Stability of an SEIR epidemic model with independentstochastic perturbations,” Physica A Statistical Mechanics and its Appli-cations , 4928–4936 (2013). L. F. Olsen, G. L. Truty, and W. M. Schaffer, “Oscillations and chaos inepidemics: a nonlinear dynamic study of six childhood diseases in copen-hagen, denmark,” Theor Popul Biol , 344–370 (1988). napshot attractors in SIR 9 D. A. Rand and H. B. Wilson, “Chaotic stochasticity: a ubiquitous sourceof unpredictability in epidemics,” Proceedings of the Royal Society B ,179–184 (1991). B. M. Bolker and B. T. Grenfell, “Chaos and Biological Complexity inMeasles Dynamics,” Proceedings of the Royal Society of London Series B , 75–81 (1993). M. J. Keeling and B. T. Grenfell, “Disease extinction and communitysize: Modeling the persistence of measles,” Science , 65–67 (1997),https://science.sciencemag.org/content/275/5296/65.full.pdf. W. P. London and J. A. Yorke, “Recurrent outbreaks of measles, chicken-pox and mumps. i. seasonal variation in contact rates,” American JournalEpidemiology , 453–468 (1973). J. L. Aron and I. B. Schwartz, “Seasonality and period-doubling bifurca-tions in an epidemic model,” Journal of Theoretical Biology , 665 – 679(1984). I. B. Schwartz, “Multiple stable recurrent outbreaks and predictability inseasonally forced nonlinear epidemic models,” J. Math. Biology , 347–361 (1985). R. M. Anderson and R. M. May,
Infectious Diseases of Humans: Dynamicsand Control (Oxford Science Publications, Oxford, 1992). S. Altizer, A. Dobson, P. Hosseini, P. Hudson, M. Pascual, and P. Rohani,“Seasonality and the Dynamics of Infectious Diseases,” Ecology Letters ,467–84 (2006). L. Stone, R. Olinky, and A. Huppert, “Seasonal dynamics of recurrentepidemics,” Nature , 533–536 (2007). G. Rozhnova and A. Nunes, “Fluctuations and oscillations in a simple epi-demic model,” Physical Review E , 041922 (2009), arXiv:0810.4683 [q-bio.PE]. S. M. O’Regan, T. C. Kelly, A. Korobeinikov, M. J. O’Callaghan, A. V.Pokrovskii, and D. Rachinskii, “Chaos in a seasonally perturbed SIRmodel: avian influenza in a seabird colony as a paradigm,” J Math Biol. , 293–327 (2013). P. G. Barrientos, J. A. Rodriguez, and A. Ruiz-Herrera, “Chaotic dynam-ics in the seasonally forced sir epidemic model,” Journal of mathematicalbiology , 1655–1668 (2017). J. Duarte, C. Januaário, N. Martins, S. Rogovchenko, and Y. Rogovchenko,“Chaos analysis and explicit series solutions to the seasonally forced sirepidemic model,” Journal of Mathematical Biology , 2235–2258 (2019). S. Bilal, B. K. Singh, A. Prasad, and E. Michael, “Effects of quasiperiodicforcing in epidemic models,” Chaos , 093115 (2016). I. Papst and D. J. D. Earn, “Invariant predictions of epidemic patterns fromradically different forms of seasonal forcing,” Journal of the Royal SocietyInterface (2019). A. Hastings and K. Higgins, “Persistence of Transients in Spatially Struc-tured Ecological Models,” Science , 1133–1136 (1994). A. Hastings, “Transients: the key to long-term ecological understanding?”Trends in Ecology and Evolution , 39–45 (2004). B. K. Singh, P. E. Parham, and C.-K. Hu, “Structural Perturbations to Pop-ulation Skeletons: Transient Dynamics, Coexistence of Attractors and theRarity of Chaos,” PLoS ONE , e24200 (2011). A. Morozov, K. Abbott, K. Cuddington, T. Francis, G. Gellner, A. Hastings,Y.-C. Lai, S. Petrovskii, K. Scranton, and M. L. Zeeman, “Long transientsin ecology: Theory and applications,” Physics of Life Reviews , 1–40(2020). D. J. D. Earn, P. Rohani, B. M. Bolker, and B. T. Grenfell, “A SimpleModel for Complex Dynamical Transitions in Epidemics,” Science ,667–670 (2000). K. Hempel and D. J. D. Earn, “A century of transitions in new york city’smeasles dynamics,” Journal of The Royal Society Interface , 20150024(2015). J. P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange attrac-tors,” Reviews of Modern Physics , 617–656 (1985). E. Ott,
Chaos in dynamical systems , 2nd ed. (Cambridge University Press,Cambridge, 2002). M. Gruiz and T. Tél,
Chaotic dynamics , 1st ed. (Cambridge UniversityPress, Cambridge, 2006). For completeness, it should be mentioned that a generalization of snapshotattractors has been done and is referred to as pullback attractors in mathe-matical and climate research. . F. J. Romeiras, C. Grebogi, and E. Ott, “Multifractal properties of snapshot attractors of random maps,” Physical Review A , 784–799 (1990). L. Arnold,
Random dynamical systems , 1st ed. (Springer-Verlag, Berlin,Heidelberg, 1998). Y.-C. Lai, U. Feudel, and C. Grebogi, “Scaling behavior of transition tochaos in quasiperiodically driven dynamical systems,” Physical Review E , 6070–6073 (1996). J. Jacobs, E. Ott, T. Antonsen, and J. Yorke, “Modeling fractal entrainmentsets of tracers advected by chaotic temporally irregular fluid flows usingrandom maps,” Physica D Nonlinear Phenomena , 1–17 (1997). Z. Neufeld and T. Tél, “Advection in chaotically time-dependent openflows,” Phys. Rev. E , 2832–2842 (1998). G. Károlyi, T. Tél, A. P. de Moura, and C. Grebogi, “Reactive Particles inRandom Flows,” Phys. Rev. Lett. , 174101 (2004). J. C. Sommerer and E. Ott, “Particles floating on a moving fluid - A dynam-ically comprehensible physical fractal,” Science , 335–339 (1993). M. Vincze, I. D. Borcia, and U. Harlander, “Temperature fluctuations ina changing climate: an ensemble-based experimental approach,” ScientificReports , 254 (2017), arXiv:1702.07048 [physics.flu-dyn]. G. Drótos, T. Bódai, and T. Tél, “Quantifying nonergodicity in nonau-tonomous dissipative dynamical systems: An application to climatechange,” Phys. Rev. E , 022214 (2016). T. Bódai and T. Tél, “Annual variability in a conceptual climate model:Snapshot attractors, hysteresis in extreme events, and climate sensitivity,”Chaos , 023110 (2012). J. D. Daron and D. A. Stainforth, “On predicting climate under climatechange,” Environmental Research Letters , 034021 (2013). G. Drótos, T. Bódai, and T. Tél, “Probabilistic Concepts in a ChangingClimate: A Snapshot Attractor Picture*,” Journal of Climate , 3275–3288(2015). M. Herein, J. Márfy, G. Drótos, and T. Tél, “Probabilistic Concepts inIntermediate-Complexity Climate Models: A Snapshot Attractor Picture,”Journal of Climate , 259–272 (2016). S. Pierini, M. Ghil, and M. D. Chekroun, “Exploring the Pullback Attrac-tors of a Low-Order Quasigeostrophic Ocean Model: The DeterministicCase,” Journal of Climate , 4185–4202 (2016). T. Tél, T. Bódai, G. Drótos, T. Haszpra, M. Herein, B. Kaszás, andM. Vincze, “The Theory of Parallel Climate Realizations,” Journal of Sta-tistical Physics (2019), 10.1007/s10955-019-02445-7. M. J. Keeling and T. D. Eames, “Networks and epidemic models,” Journalof the Royal Society, Interface , 295–307 (2017). R. M. Anderson and R. M. May, “Directly Transmitted Infectious Diseases:Control by Vaccination,” Science , 1053–1060 (1982). D. L. In Epidemiology: SIMS 1974 Utah Conference Proceedings andS. K.L. Cooke (eds.), eds.,
Transmission and control of arbovirus diseases (Philadelphia, 1975). G. Katriel and L. Stone, “Attack rates of seasonal epidemics,” Math Bio-science , 56–65 (2012). E. N. Lorenz, “Irregularity: a fundamental property of the atmosphere,”Tellus Series A , 98 (1984). Choosing a different initial phase does not effect the qualitative picture ofthe dynamics. Only the shape of the phase space object differ to those ob-tained in present study. A. Becker, A. Wesolowski, O. Bjørnstad, and B. Grenfell, “Long-termdynamics of measles in london: Titrating the impact of wars, the 1918 pan-demic, and vaccination,” PLoS Comput Biol , e1007305 (2019). B. Kaszás, U. Feudel, and T. Tél, “Death and revival of chaos,” Phys. Rev.E , 062221 (2016). Y.-C. Lai and T. Tél,
Transient chaos /Lai-Tel. Berlin : Springer,c2011. 2011 (Springer, 2011). Generally speaking, the time variability of b always results in a changingphase space structure according to the theory of snapshot attractors. M. Ghil, M. D. Chekroun, and E. Simonnet, “Climate dynamics and fluidmechanics: Natural variability and related uncertainties,” Physica D Non-linear Phenomena , 2111–2126 (2008), arXiv:1006.2864 [math.DS]. M. D. Chekroun, E. Simonnet, and M. Ghil, “Stochastic climate dynam-ics: Random attractors and time-dependent invariant measures,” Physica DNonlinear Phenomena , 1685–1700 (2011). napshot attractors in SIR 10