Standard epidemic models based on compartmental differential equations are investigated under continuous parameter change as external forcing. We show that seasonal modulation of the contact parameter superimposed a monotonic decay needs a different description than that of the standard chaotic dynamics. The concept of snapshot attractors and their natural probability distribution has been adopted from the field of the latest climate-change-research to show the importance of transient effect and ensemble interpretation of disease spread. After presenting the extended bifurcation diagram of measles, the temporal change of the phase space structure is investigated. By defining statistical measures over the ensemble, we can interpret the internal variability of the epidemic as the onset of complex dynamics even for those values of contact parameter where regular behavior is expected. We argue that anomalous outbreaks of infectious class cannot die out until transient chaos is presented for various parameters. More important, that this fact becomes visible by using of ensemble approach rather than single trajectory representation. These findings are applicable 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
The effort to estimate short and long-term behavior, system parameters, possible control strategies, and also social impact of COVID-19 pandemic stimulated an explosion of disease modelling 1 and scientifically demanding examinations from many different point of view made by various community members and research groups around the globe. Although the present situation is rather new, we have many examples from the history that helped to lay down basic mathematical rules of disease propagation. The celebrated SIR (susceptible -infectious -recovered) model 2 splits the population into three disjoint parts and deals with the number of individuals in 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 good agreement with observations. In a homogeneous environment these models possess a globally stable fixed point attractor as a disease equilibrium 3 . Omitting the stochastic nature of real world disease spread [4] [5] [6] [7] , the deterministic SIR-like models present diverse and reach dynamics. The non-linearity of the model [8] [9] [10] [11] and also the time dependent internal forcing can be considered as a source of complexity. Sufficiently large seasonal forcing or different mixing rates between sub-populations can cause large oscillations and also period doubling cascades [12] [13] [14] [15] [16] [17] [18] . It has also been demonstrated that for certain values of system parameters the nonautonomous models of recurrent epidemics (measles, mumps, rubella, H5N1 avian influenza) show chaotic behavior 6, [19] [20] [21] . Traditional SIR-like epidemic models are dissipative nonlinear low-dimensional systems with constant, periodic, quasiperiodic 22 , or term-time 3,23 external forcing whose dynamics a) Electronic mail: tkovacs@general.elte.hu is governed by (chaotic) attractors in phase space. Furthermore, even if the long-term dynamics is regular and the final state of the system is a fixed-point attractor the route to this condition might be rather complex. Many studies examined the role of finite-time irregularity in ecological models [24] [25] [26] [27] as well as in epidemic dynamics 9, 23, 28, 29 concluding the relevance of transient behavior.
It is known from dynamical systems theory that chaotic attractors always can be characterized by the natural measure corresponding to the distribution of possible states in phase space [30] [31] [32] . It turned out, however, that traditional numerical methods, such as monitoring a single chaotic long-term trajectory, fail in case of arbitrary forcing (like smoothly varying parameters). Obviously, if one wants to model such a system (e.g. the changing climate with shifting atmospheric CO 2 concentration), this shortcoming has to be overcome. The mathematical concept of snapshot attractors 33 , known for many years in (theoretical and experimental) dynamical systems community [34] [35] [36] [37] [38] [39] [40] [41] , fulfills entirely our wish.
Briefly, snapshot attractors are time-dependent objects in phase space of non-autonomous dynamical systems. Thus the shape of snapshot attractors is changing in time while their fractal dimension might even remain constant 39 . Furthermore, obtaining one single trajectory in a system with arbitrary driving force, does not provide the same result as the ensemble approach (many trajectories emerging from slightly different initial conditions) in the same system at a given time instant. This effect is the consequence of the fact that ergodicity is not satisfied as the system is driven aperiodically 42 . This conclusion generated the recent opinion that the changing climate should be scrutinized by ensemble approach (parallel climate realizations) rather than averaging a single long-term time series [43] [44] [45] [46] [47] [48] .
In this work a seasonally forced deterministic epidemic model with monotonically changing contact rate (due to, for instance, vaccination or restricting the movement of the popu-lation) is presented. We propose that the statement by Ref. 43 "climate change can be seen as the evolution of snapshot attractors" also holds for disease spread dynamics with continuously 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.
Compartmental disease models describe the number of individuals in a population regarding their disease status: susceptible (S), infectious (I), or recovered (R). Although, these models involve many simplifications (such as the progression of infection or difference in response of individuals) they performed well in real world epidemic situations 49 . There are two major groups of epidemic models: the SIR-like cluster that characterizes lifelong immunity (e.g. measles, whooping cough), and the SIS-like class (containing mostly sexually transmitted diseases) which portrays repeated infections 49 .
Here we study the SEIR equations 50 that involves a fourth group in addition to previous ones. More concretely, we assume that an individual enters the population at birth as susceptible and leaves it by death. A susceptible becomes exposed (E) when contacting one or more persons, called infective(s), that can transmit the disease. In an incubation period the exposed individuals are infectious but are not yet infectious. After this term they become infective and become later immune or recovered.
The mathematical model associated with above description reads as follows
with the following notations and assumptions 51 : (i) S, E, I and R are smooth functions of time and the size of the whole population remains unchanged, S + E + I + R = 1; (ii) there are equal birth and death rates (m); (iii) the probability that an exposed individual remains in this class for a time period τ after the first contact is exp(−aτ), where 1/a is the mean (or characteristic) latent period; (iv) similarly 1/g gives the mean infectious period, the time period that an individual spends as infectious before recovery; (v) immunity is permanent and recovered individuals do not re-enter to susceptible class. The contact rate b(t), the average number of susceptibles contacted per a single infective per unit time, is the origin of the spread of the disease. In the case of annual periodicity
and t is measured in units of years.
When the system parameters are kept constant (m, a, g, b 0 , b 1 = 0), the solution of Eq. (1) shows weakly dumping oscillation with a globally stable equilibrium (b > g). In this case the dynamics can be further characterized by the expected number of secondary cases caused by an infectious individual in the susceptible population. This is known as the basic reproductive ratio, and can be expressed denoted here by R 0 = ba/[(m + a)(m + g)]. 13 R 0 essentially determines whether a disease can (R 0 > 1, endemic equilibrium) or cannot (R 0 < 1, infection-free) persist in a population. Provided that latent and infectious periods are short, i.e. a, g m, one can use the approximation R 0 = b/g.
There are observations, for example in childhood diseases or avian influenza, that do not show damped characteristics rather (ir)regular cycles instead. This phenomenon can be linked to seasonal variation in contact rate b(t) or in recruitment rate g 19 . It has already been noticed that in case of periodic forcing, b 0 (t) = b 0 =const. in Eq. (2), chaos emerges in epidemic time series 3, 21, 52 . Similarly to the Lorentz84 model 53 wherein the solar irradiation induces seasonal effect, the SEIR model with variable contact rate is a nonautonomous low-order system with well-defined chaotic attractor in the phase space at certain parameter values. The low dimension of the phase space allows to visualize its pattern comfortably.
As a result of the fixed population size, point (i) above, one of the four equations (traditionally the equation for R) can be omitted. In addition, the infectious and exposed groups turn out to be linearly related to the first order 13 at least, see Figure S1 in supplementary material. Consequently, the S − I phase portrait represents the state space texture accurately. Plotting one trajectory in the S − I plane of the phase space one would observe a coil design due to the explicit time dependence of the model, b 1 = 0. A conventional method to make a periodically forced system autonomous is to take "pictures" about the phase space with the same frequency as the excitation acts, i.e. making a stroboscopic map 31,32 . By choosing an initial time t 0 and purely sinusoidal force in Eq. 2 (b 0 (t) =const.) one can interpret the filamentary shape of the chaotic attractor by defining the stroboscopic map t 0 mod T = 0. Therefore, the periodic variation in Eq. (2) yields am annually stationary epidemic with steady attractor.
Eq. (2) assumes a constant average contact rate b 0 . We have an annual period starting with its maximum at t = 0. 54 In physical or biological context this means the largest contact rate possibly due to school term, vacation and holidays, seasonal breeding pattern in a seabird colony, etc. Eq. (1) with the following parameters, corresponding to measles 6 , shows irregular dynamics: m = 0.02 year −1 , a = 35.84 year −1 , g = 100 year −1 , b 0 = 1800 year −1 , b 1 = 0.28 year −1 . The above parameters lead to R 0 ≈ 18. Numerical calculations show that for lower mean values of the contact rate periodic fixed-point attractors exist. For instance, biannual cycles at b 0 ≈ 1500. For more detailed asymptotic dynamics see the bifurcation diagram in Figure 1 (blue curve).
Until now b 0 was kept as a constant. Following the climate change methodology monotonic variation in b(t) through b 0 (t) is established. Equation (3) specifies the temporal de-cay of the mean contact rate
Here t st represents the time until the epidemic is stationary, i.e. b 0 = const., this is set to be 250 years after initialization. This amount of time seems to be enough since the convergence of trajectories to the attractors is much shorter, t c ≈ 50 years. The exponent α is the decay rate of b 0 (t). In this study α takes three values 0,04 yr −1 , 0.01 yr −1 , and 0.004 yr −1 , and b 0 always starts at 1800, i.e. from the chaotic regime. Eq. (3) implies that lim t→∞ b 0 (t) = 90. This means that the critical value of mean contact rate b 0 = 100 (or in terms of basic reproduction rate R 0 ≈ 1) is reached at different times after t st , according to α.
The main motivation of this study is based on Eq. (3). Considering, for example, the time series of childhood measles in London after World-War II (between 1945 and 1990), the vaccination program (starting in 1968) changed dramatically the number of cases 3, 28, 55 . Beside the medical treatment, if available, other artificially forced decline in contact rate (such as gradual lockdown prescribed by administration) can be imitated by Eq. (3) resulting in a changing epidemic. This secular variation of parameter b 0 gives a clear analogy to climate change models and the framework applied in this field.
The bifurcation diagram basically reveals the long-term dynamics, excluded the initial transients, in a given parameter range. Classically, one single trajectory sampled by stroboscopic map serves the blue shape of bifurcation diagram in Figure 1 . In case of the stationary or stationary epidemic (b 0 = const.) complex dynamics arises for certain parameter values in the SEIR model, b 0 1770. We should note, and it is going to be essential in our analysis, that the bifurcation diagram can also be established in a different way. That is, a large number of initial conditions are placed in the phase space and their evolution is monitored for sufficiently long time, but much shorter (one or two orders of magnitude) than that of the single trajectory considered above. Due to the ergodicity of the stationary chaotic dynamics the end-points of the ensemble members portray exactly the same pattern 56 .
The same applies for the phase space pattern, too. Fig. 2 depicts the chaotic attractor obtained after T = 2500yr integration and 1600 initial conditions distributed uniformly in phase space. The same picture is obtained by one single trajectory after 200k year simulation. A representative S − I phase portrait is pictured in Figure 2 . The chaotic attractor (b 0 = 1800) has been established fairly early and it remains unchanged when viewed after integer multiples of T = 1yr. Nevertheless, choosing a different day of the year, the same applies eith a different pattern in S − I plane. The variable S is plotted on the same year of the day starting at t 0 = 0. The classical diagram (blue curve) is obtained by integrating 2000 initial conditions for t = 15000yrs and the last 500 points stored. For low parameter values cycles are determined by annual periodicity. There are many period doubling visible (shorter blue segments), however, the first one that routes to chaos starts at b 0 ≈ 1250. The end-points of the same trajectories are plotted after 500yr iteration indicating long-lived transients (green). Clearly finite time chaotic motion unfolds the extra structure of bifurcation diagram. Upper x-axis evaluates the basic reproductive ratio derived from system parameters. R 0 = 1 is also marked (red arrow) to indicate the edge of endemic equilibrium. Long-term dynamics for R 0 < 1 (b 0 < 100) implies a dying out infection in population. The inset display real permanent chaos. Top box displays the R 0 of common diseases for comparison only. (Ebola * : based on 2014 Ebola outbreak.) Transient chaos, complex behavior on finite-time scales 57 , manifests before a trajectory comes to an attractor, both in dissipative as well as in conservative systems. In general, the attractor can be a simple object in the phase space, for instance, a fixed-point or a limit cycle. Following the time evolution of the ensemble before it approaches the attractor, one can capture other ingredients of the bifurcation diagram by considering the initial transients. To do this, the uniformly distributed 2000 trajectories are integrated for 500 years and the corresponding state space positions are stored. This feature becomes visible in the green vague domain that extends along the blue curve. Clearly, transient chaos has a significant contribution to the dynamics for b 0 125.
The initial transients are usually thrown away in long-term dynamical analysis. However, the complexity might appear just in this phase of motion as indicated by the bifurcation diagram in Figure 1 Introducing classical statistical measures on the ensemble one can quantitatively keep track of transient effect. The mean
denotes an observable corresponding to the ith member at time t) designates how fast the initial irregularity terminates and also adverts the average transient time. One can observe the similar characteristics for both the S and I curves, fixed-point and chaotic attractors (top and middle panels), respectively.
In the stationary epidemic model one might expect that after some time the extent of the attractor remains constant suggesting that all of the individual trajectories in the ensemble have reached it. The standard deviation,
refers to the size of the attractor in the direction of A. Stationary epidemic with constant σ A essentially means that after the transients the phase space structure does not change in time. This can be visualized by using of ensemble approach. Without loss of generality we always start our simulation at t = 0 (that corresponds to that part of the year with highest seasonal amplitude b(t)) and take the next picture about the ensemble at t mod T = 0 that coincides the same day of the year.
The natural distribution of a chaotic attractor 31 comes from the fact that the dots visit certain parts of the filamentary structure more frequently than others. Moreover, this distribution is stationary and does not depend on initial conditions. Instead of investigating a 3D frequency diagram (see supplementary Figure S2 ), we plot the projection I of trajectories wandering on the chaotic attractor in Figure 2 . To obtain Figure 4 , a fine grid is defined in the S − I phase plane and the number of points in each cell is plotted against the variable I. The histogram consists of four different ensembles. Ensemble 1,2, FIG. 4. Natural distribution on the chaotic attractor (b 0 = 1800) projected onto the I variable after t = 8500 years. The distribution does not depend on the individual members of the ensemble. The frequency has been cut at 3500 for the clarity. and 3 covers the same volume in the phase space but involves slightly different initial state vectors. Ensemble 4 lays out initial conditions from other part of phase space. The distribution clearly shows the same pattern for all four ensembles.
In a changing epidemic the parameter(s) of the system is(are) varying as time goes on. One might think that, say, a decreasing of the contact parameter b 0 means to walk along the bifurcation diagram slowly from right to left, and terminating at a safe destination with R 0 < 1. In this section we will show that this idea is fairly naive due to the internal variability and the transient effects in the epidemic. To capture the dynamics properly in this scenario, the ensemble approach and the concept of snapshot attractor is desirable.
According to the common measles parameters we start the switch-off process, called also epidemic change, from the chaotic attractor (b 0 = 1800) according to Eq. (3).
Recent studies 56 reveal that starting from the chaotic attractor well after that the trajectories reached it, the evolution of different states might be rather diverse. That is, transient dynamics becomes important again while the parameter change appears in the system. In Figure 5 , four trajectories are selected from the chaotic attractor and their time dependence is followed. It can be immediately seen that different colors reach the fixed-point attractor at different times. For instance, in case of α = 0.004 (bottom curve), first the blue, then the yellow, red, and green curves arrive at the fixed-point attractor. One can pick up other trajectories that will have longer or shorter oscillations. Furthermore, the previous order of colors, i.e. transient times, might change with α, as demonstrated by the middle curve. Therefore, we can point out, that it is worth investigating several trajectories simultaneously instead of following individuals. However, after the parameter shift is switched on, the shape of the chaotic attractor starts to vary and forms a time-dependent set called snapshot attractor. Thus, snapshot attractor is an object in the phase plane that contains the whole ensemble at a given time instant. We note that this mutation is not a result of the particular choice of initial phase (day of the year) of the mapping rather than the change of parameter. The geometrical alteration of the snapshot attractor is followed by the change of distribution on it. Nevertheless, its form and the distribution is independent of the choice of the ensemble and the onset (t st ) of parameter decay. Figure 6 exhibits snapshot attractors drawn by N = 10 4 initial conditions during the epidemic change. The mean contact rate decreases from top to bottom and left to right as specified on each panel. Thus, every single plot corresponds to a different time instant of simulation. The top left panel, taken at T = 219yr, coincides with the attractor in Figure 2 since the parameters of the system are the same for t < t st = 250 yr. A more interesting aspect shows up on the following two panels. A filamentary structure, the fingerprint of chaotic behavior, still dominates the pattern although the parameter b 0 is well below of its original value (1800) as well as of 1770 where large extended chaotic attractors are formed in the bifurcation diagram (Fig. 1 inset) . In other words, for those contact rates (b 0 = 835, 328) no chaotic behavior is anticipated in the stationary dynamics.
Although, Figure 6 indicates the transmutation of the snapshot attractor for decay exponent α = 0.01 one can obtain similar alteration for different switch-off rates too. Slower scenario (e.g. α = 0.004 in Eq. (3)) allows the attractor to keep its filamentary shape and maintain chaotic dynamics longer. From other perspective, the same parameter value b 0 is reached later while α is smaller. The opposite is true for a faster parameter change, say α = 0.04. It can also be shown that approximately the same pattern belongs to the same contact value regardless the rate of change. That is, the faster the contact rate decay, the less pronounced the complexity in phase space. Other interesting feature is that although the contact rate b 0 decreases monotonically, the size of the snapshot attractors might increase, Figure 6 top right panel. This fact yields that the domain accessible by the dynamics can be larger. Thus the possible (S, I) pairs extend to larger domain in phase space even for smaller contact parameter. Furthermore, the average (green circles) and the standard deviation (red ellipses) also change in time. And this temporal behaviour cannot obtained from the classical view, only by using snapshot attractors.
To understand this phenomenon we should recall the concept of chaotic saddle. This non-attracting set with its stable and unstable manifolds in phase space is responsible for the finite time chaotic behavior 57 . To construct the saddle numerically, we define two holes in the phase space (black rectangles in Figure 7) . Then, a large number of initial conditions (N = 3 · 10 5 ) are distributed uniformly in the region S = [0.04; 0.2], E = [0; 0.00005], I = [0.00001; 0.01] and the trajectories a integrated for T = 30 years. Only those trajectories are kept that do not enter the holes during the integration time. The initial conditions belonging to these trajectories draw the stable manifold of the saddle (red dots in Fig. 7) , while those that are just before leave plot the unstable manifold (green dots.) The saddle itself is the intersection of its manifolds, not shown here. The average lifetime of chaos (the inverse of escape rate, κ) can be estimated by calculating the time distribution of the non-escaped trajectories.
It is also known from the theory of transient chaos that the saddle's manifolds have filamentary design just like the snapshot attractors in Fig. 6 . Analogously to chaotic attractor, stationary dynamics with constant driving amplitude defines a stationary chaotic saddle related to certain parameter values of the system. Due to the continuous adjustment of the contact rate, b(t), in the changing epidemic model, the trajectrories do not have time to reach the attractor belonging to a given b 0 value. Consequently, a time dependent chaotic saddle is considered, whose unstable manifold approximates the snapshot attractor 56 . This stationary saddle persist for very low b 0 values and its unstable manifold controls the finite time complex epidemic dynamics, Figure 7 .
We emphasize at this point that the natural distribution associated to the snapshot attractor also changes in time following the geometrical reorganization of the phase space pattern (histogram visualization in Fig. 6 ). Similarly to climate research we can define the concept of parallel epidemic realizations. In standard disease models like Eq. (1) this can be done by using a large number of initial conditions as an ensemble in phase space and following their evolution as discussed in Section III A. This picture can be imagined as many copies of the epidemics obeying the same physical laws and being affected by the same time-dependent forcing 48 .
As we have seen before, in case of a changing epidemic, the time-dependent forcing has an impact on the natural measure of the snapshot attractors. This fact results in a temporal change of the average values as well as internal variability. To quantify the internal variability of the system (1) statistical measures over the ensemble should be proposed like in the case of stationary epidemic. The variance or the standard deviation σ A in Eq. (4) of the ensemble characterize the fluctuations around the averages indicating the inherent internal variability.
Top panels of Figure 8 show the ensemble standard deviation σ S (blue) and σ I for different decay rates α. The parallel epidemic realizations contain 1600 numerically integrated trajectories sampled at integer multiplies of years. The convergence time t c , the time needed of the ensemble to reach and spread along the attractor, is found to be t c = 50yr. Both σ S and σ I are constant before t st = 250yr since till then the seasonal driving is constant, i.e. b 0 = 1800yr −1 . Right after this the contact rate starts to decrease (t > 250yr) according to Eq. (3) both the mean state and the internal variability of the epidemic changes with time.
In changing epidemic, first, the graph of σ S (t) increases, after the maximum the trend follows nearly b 0 (t), it decays to zero illustrating that the size of the attractor shrinks and asymptotically reaches the neighbourhood of a regular fixedpoint attractor, (S, I) = (1, 0), as expected for b 0 ≈ 90 (R 0 = 0.9). The larger the α, the more regular, i.e less filamentary, the phase space geometry at the same time, see bottom panels in Fig. 8 . The shape of the standard deviation can be explained by transient chaos that occurs during the parameter change. The size of the green shaded band around the main feature in the bifurcation diagram (Fig. 1) already indicates that the size of the phase space region filled with transient chaos increases as b 0 reduces. To demonstrate this we call the attention to the horizontal dimension of the snapshot attractors in Figure 6 . The smoothly decreasing profile of the standard deviation after reaching the maximum refers to the existence of transient effect up to very small parameter values.
Different maximum values of the σ S (t) curves indicate a non-trivial relation between the internal variability and changing rate of b 0 . The largest value corresponds to α = 0.01 (top middle panel), the other two scenarios (α = 0.04 and 0.004) show roughly the same maximum size of the snapshot attractor, albeit, at different time instants. Bottom panels present the maximal physical extension of the snapshot attractors versus time. These plots also supports our observations that the size of the attractor first increases and the shrinks to be a fixedpoint attractor.
One possible explanation of this property is the relative "coupling" between the timescales, that is, decay of forcing defined by the exponent α and the average lifetime of chaos in stationary dynamics at specific parameter b 0 . From a geometrical point of view, how close the snapshot attractor evolves to the unstable manifold of the stationary chaotic saddle at particular contact parameter. A detailed exploration of this feature is, however, postponed to a future study.
Parallel climate realizations, an effective and new framework in climate research, has been adopted to an epidemic model to explore the fading of the complexity due to the systematic switch off the driving mechanism. The mathematical concept of snapshot attractors and their natural distribution demonstrates that single time series analysis is not capable to reflect the complex dynamics of a changing epidemic. Instead of monitoring isolated events, the ensemble view of trajectories -parallel epidemic realizations -and its statistical description is desirable.
The temporal change of the attractor geometry and the distribution on it reveals the internal variability of the dynamics.
The extension of the bifurcation diagram indicates the importance of transient chaos in the stationary dynamics as well as during the continuous parameter shift. No matter whether we start from a stationary state of the long-term dynamics, the switch-off process activates the hidden parts of the bifurcation diagram. Thus, the underlying non-attracting object, the chaotic saddle, or more precisely its unstable manifold, orga- nizes the system's evolution.
The dissipative relaxation timescale, the inverse of the phase space contraction rate based on the divergence of the vector field of system (1), is extremely low, ∼ 2 − 5 days, for measles. The other characteristic times are, the inverse of the switch-off rate and the escape rate from the chaotic saddle are, for comparison, α −1 ≈ 25 − 250 yrs, κ −1 ≈ 30 − 35 yrs, respectively. 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 56 . In other words, there is not enough time for trajectories to reach the stationary attractor, either it is chaotic or regular, due to the nonstop parameter variation. Therefore, the dynamics of the switching off process shows much more complexity than that of the stationary scenario.
In the spirit of the investigated SIR-like classical epidemic model (1) with parameter shift 58 of the contact rate b(t) we emphasize the importance of the finite time irregular dynamics. Our results reinforce the "hidden" complex transients, and also the importance of time-dependent snapshot attractors not just for an epidemic with large reproductive ratio, such as measles and chickenpox, but for those with lower values (R 0 1) too.
The main conclusion of this study sheds light on the importance of the ensemble view and parallel realizations in epidemic dynamics. Fig. 6 illustrates the general features one expects in the (long-term) prediction of any epidemics. Predictions based on individual simulations are not reliable since, due to the chaotic nature of the dynamics, they can lead to many possible results, to any point of the snapshot attractor belonging to the time instant of the prediction. The ensemble approach is able to treat all possible outcomes as a whole, and predicts even the probability of the different permitted epi-demic states. When convenient, statistical moments of this distribution can be determined, e.g. the average (providing the most typical epidemic outcome), or the variance (characterizing how broad the distribution of the permitted states is), but higher order moments might also be useful. Anyhow, one can thus follow how the statistical prediction changes in time.
We benefited from useful discussions with T. Tél. This work was supported by the NKFIH Hungarian Grants K125171. The support of Bolyai Research Fellowship and ÚNKP-19-4 New National Excellence Program of Ministry for Innovation and Technology is also acknowledged.
FIG. S1. Three dimensional chaotic attractor in the stationary SEIR model. The parameters are defined in the main text. It is striking that the infectious (I) and exposed (E) classes are approximately linearly related 13 (according to the ratio g/a). The relation is also indicated with red dashed line in the I − E surface. Thus, one can investigate the S − I projection without loss of generality. This has been done in the present work.
FIG. S2. Natural distribution on the chaotic attractor of Fig. S1 . A fine grid is defined in the S − I plane and the number of phase points in each cell is recorded. The histogram shows that some part of the attractor are more frequently visited than others assigning the natural measure on it. Note that the height of the bars along the I ≈ 0 axis is 1-2 order of magnitude larger than in the rest of the plot.
We present two illustrative plots related to the phase space structure of SEIR model (1) that complement the overall picture of the paper but lack of the plots in the main text does not violates the understanding of the basic idea.