The surveillance data we worked with were the weekly shares of positive influenza tests for all influenza types in the Northern Hemisphere from January 5, 2009, to March 31, 2025.

The data source is from Our World in Data, which processed the data from FluNet by the World Health Organization (2023), and the data are available under a Creative Commons license at the Our World in Data website [26].

The data are an aggregate indicator on the share of positive tests for the Northern Hemisphere calculated by the Our World in Data team from the FluNet original source and are part of the influenza data tracker made available at the Our World in Data website for researchers.

Even though the actual number of potential cases may individually be underestimated and country-specific data vary in accordance with the testing frequencies and protocols, the aggregate number calculated by the Our World in Data for the Northern Hemisphere as a whole mitigates, to some extent, the different testing frequencies and protocols by different countries. This makes it a surveillance series that is correlated with actual influenza outbreak processes and can provide a useful first-start, exploratory surveillance data series for studying the possibility of complex chaotic and stochastic epidemiological dynamics in influenza. This allows the application of the predictive decomposition method with epidemiological interpretability, as we show in both the Results and Discussion sections.

The analysis followed a 3-stage workflow (Figure 1), which we describe in the following paragraphs.

The predictive phase space reconstruction method using topological machine learning is a critical step and is explained in detail in Multimedia Appendix 1. In this method, the observed time series was transformed into a geometric representation of system states using delay embedding in a way that maximizes prediction [27].

The reconstructed trajectory was analyzed to determine whether it exhibited evidence of chaos or random behavior. This step established whether the observed dynamics contained meaningful structure or were primarily driven by noise.

The series was decomposed into a predictable component and a stochastic component. This step was necessary because epidemiological dynamics combine structured dynamics with random fluctuations. The decomposition was carried out until the residuals behaved as IID noise, meaning that no systematic linear nor nonlinear patterns remained, and the remaining variability was effectively random.

The weekly share of positive influenza tests showed an intermittent pattern, with long low-activity periods followed by sharp outbreak peaks (Figure 2). This indicates that influenza dynamics are not uniform over time but instead occur in bursts consistent with epidemic waves.

The power spectrum exhibited a power-law decay with a slope of −2.52 (R²=69.3%, P<.001; Figure 3). This showed that the series has strong persistence, meaning past values influence future behavior over extended periods [11,28]. This last result excluded the hypothesis of white chaos and indicated that, if chaos is present, it will be a form of color chaos [16].

Using a k-nearest neighbors’ regressor with a 10-week rolling window, the optimal embedding was obtained for dimension 2 and lag 11, achieving R²=78.33%, explained variance of 78.34%, correlation with the target series of 0.8902, and a root mean square error (RMSE) of 4.7277, which represented 11.23% of the series’ amplitude (Table 1). This indicated that a low-dimensional representation is sufficient to capture most of the system’s dynamics in a way that can support epidemiological prediction.

To test the robustness of underlying 2-dimensional deterministic dynamics linked to influenza that can be used to predict the target series with a high prediction score, we needed to plot the optimal dimensions obtained for different lags (Figure 4) and, for the optimal lag, plot the score for different dimensions (Figure 5).

From Figure 4, we concluded that prediction performance exceeded 77% across all tested embeddings, which showed strong predictability using phase space embedding. We also found that dimension 2 was the most frequently selected, which demonstrated the robustness of the 2-dimensional phase space representation.

Performance decreased as embedding dimension increased beyond 2 to 3 (Figure 5); therefore, increasing the embedding dimension did not improve representation and may have reduced predictive accuracy. This indicated that the 2-dimensional structure was better at capturing the underlying epidemiological dynamics.

Figure 6 shows the phase portrait of the trajectories. The reconstructed attractor exhibited a conic shape with a concentration of points near the origin. This showed that the system spends extended periods in low-activity states, with more dispersed trajectories corresponding with higher-activity phases and with the epidemiological outbreaks.

In Figure 6, there is an overall coincidence of the original reconstructed attractor and the k-nearest neighbors’ predicted trajectory, which further reinforced the aforementioned findings of strong predictability and that there was a strong deterministic feature to the dynamics.

The largest Lyapunov exponent (LLE), estimated using the method by Eckmann et al [29], was positive for both the original and predicted attractors (7.1^–04 and 9.1^–04, respectively). The Kaplan-Yorke dimension lies between 1 and 2 (1.8587 for the original attractor and 1.9778 for the predicted attractor). The box-counting dimension estimates were in close agreement with the Kaplan-Yorke dimension (1.7891 and 1.8032, respectively; P<.001), which confirmed the presence of a fractal chaotic attractor.

Both estimated LLEs were statistically significant at the 1% level based on surrogate testing using the reconstructed, predicted, and wavelet-filtered attractors for the power law surrogate, excluding a power law noise process. For the stochastic version of the Iterative Amplitude Adjusted Fourier Transform (SIAAFT), the original reconstructed attractor’s LLE was statistically significant at the 5% level, with both filtered attractor’s LLEs statistically significant at the 1% level (Table 2).

These results strongly reinforced the evidence that we were dealing with a chaotic attractor, rather than a purely random process. This implied that the dynamics retained an underlying deterministic structure that can be used for short- to medium-term prediction.

Using the rolling window k-nearest neighbors’ predictions as the benchmark, we then compared the performance of different major machine learning–based predictive algorithms on a 10-day rolling window [42,43].

Ridge regression (λ=0.5) achieved the highest performance for 1-week-ahead prediction (R²=92.11%, RMSE=2.85; Table 3). This showed that short-term dynamics can be accurately approximated using rolling window linear models (Table 3).

For multistep prediction, gradient boosting performed best for longer horizons (Table 4). This indicated that more complex nonlinear patterns became relevant as the prediction horizon increased. Prediction accuracy remained greater than 73% at 5 weeks (Table 4). This showed that the system retained predictive structure beyond short-term horizons.

Residuals from the ridge regression model showed clustering of large fluctuations in bursts (Figure 7), indicating time-varying volatility. The Engle test confirmed the presence of autoregressive conditional heteroskedasticity (ARCH) effects (Table 5), meaning that periods of high and low variability clustered over time rather than occurring randomly.

After applying an exponentially weighted moving average model (λ=0.1; Multimedia Appendix 2), no ARCH effects remained (Table 5). This indicated that volatility was largely captured by the model and occurred in short-lived bursts rather than persistent regimes [35-40].

The standardized residuals passed the BDS test for independence (Table 6) [41]. This showed that the remaining noise behaved as an approximately independent process.

The residual distribution was consistent with a Student t distribution (Table 7). This indicated the presence of heavy tails, meaning extreme deviations were more likely than under a normal distribution.

None of the major generalized ARCH models estimated on the full sample passed the BDS test (Table S1 in Multimedia Appendix 2), which indicated that the exponentially weighted moving average model performed better at capturing the volatility pattern.

The estimated volatility exhibited burst-like increases over time (Figure 8A). This indicated that uncertainty rises during periods of higher influenza activity.

Outbreak probability varied over time as a function of both predicted levels and volatility as exemplified in Figure 8B for an outbreak event probability associated with a share of positive tests higher than 15. This showed that risk depended not only on expected activity but also on the level of uncertainty.

Periods of high influenza activity coincided with higher volatility. This suggests that outbreaks are associated with increased dispersion and elevated epidemiological risk.

In this work, we were able to decompose the dynamics of the weekly share of positive tests for influenza in the Northern Hemisphere from January 1, 2009, to March 31, 2025, into a component that is driven by a low-dimensional chaotic attractor characterized by black chaos with long-range persistence and a stochastic component that is the result of a time-varying volatility parameter multiplied by IID noise.

We now discuss the implications of the main results on 3 levels: the first level regards the connection to previous literature on theoretical epidemiological modeling and the implication in terms of reinforcing this literature and providing new directions for theoretical epidemiological modeling; the second level regards the connection to the previous empirical work on chaos in epidemiology including in influenza surveillance data; and the third level regards the implications for epidemiological risk analysis, monitoring, and surveillance.

Considering the first level, our results reinforced the theoretical studies that have found chaotic dynamics to be possible in basic epidemiological models [1-8], while these models raised the hypothesis of chaos in general epidemiological dynamics, including the issue of whether chaos can actually be empirically present in these dynamics and what types of chaos are a key problem for empirical research on epidemiology.

A basic implication of chaos for epidemiology is that policymaking may have nontrivial, nonlinear consequences on the dynamics, including the possibility of triggering chaotic dynamics and turbulent bursts as addressed in [9], or even reduce them, as took place during the COVID-19 lockdown for the influenza series that we analyzed.

Another implication of chaos is that it has a complex relationship with predictability, in the sense that, although chaos implies limits to predictability due to the sensitive dependence upon initial conditions, chaotic attractors also have topological signatures that allow for predictability. In addition, there are different types of chaos, ranging from white chaos, which has a white noise–like spectrum and is less predictable, to color chaos, which has long-range dependence that makes it more predictable.

In our empirical results, the evidence was favorable to the hypothesis of color chaos in the form of black chaos, which is consistent with epidemiological models that exhibit chaotic intermittency associated with epidemiological seasonality and with previous studies using influenza data [19].

Although consistent with the theoretical possibility of chaos in epidemiological models, the decomposition down to IID noise surfaced an issue regarding these models. A major mathematical approach to the modeling of stochastic chaos usually involves adding some form of IID Gaussian or uniform noise to a deterministic set of equations that can exhibit deterministic chaos [20,22]. The findings from our decomposition method, however, did not support this assumption in the epidemiological context.

Indeed, we found that the stochastic dynamics associated with the noise factor in the predictive decomposition, rather than IID, showed a nonstationary variance with ARCH-like dependence. This means that epidemiological mathematical modeling that produces chaos in SIR and SEIR models needs to consider both the impact of simple IID noise. It also needs to compare it with what happens when stochastic processes with ARCH dependence in the variance are introduced, since these processes may lead to significant qualitative changes in the chaotic dynamics, namely, in the connection between outbreaks and turbulence.

Through the decomposition method, we found that outbreak waves visible in the surveillance series under study tended to occur in tandem with a rise in variance. This leads to a rise in dispersion in these periods as well as turbulence bursts associated with the outbreak waves.

The other impact of the ARCH stochastic component is that it injects short-run fast turbulence in the attractor’s long-range dynamics, leading to stochastic turbulence leaks into the chaotic dynamics.

The study of these phenomena in stochastic versions of general theoretical epidemiological models, like SIR and SEIR models, that exhibit deterministic chaos is a relevant research direction that may feed back into new empirical methods for epidemiological simulation, policymaking, and risk modeling.

The second dimension of this work regards its connection to previous empirical works on chaos in epidemiological surveillance data. We found evidence of a low-dimensional fractal attractor with a positive LLE that was low but statistically significant upon surrogate testing on both the original and the noise-filtered series. The exponent after wavelet noise filtering was close to those found for COVID-19 surveillance data, though slightly smaller [11-14].

Regarding the influenza data, our work reinforced the findings in [19], in which statistically significant (upon surrogate testing) evidence of chaos was found in influenza data, with positive LLEs and long memory consistent with color chaos, which is also consistent with our findings.

However, by decomposing the attractor, we also found additional evidence regarding the dynamical noise affecting the influenza epidemiological dynamics that is nonstationary in variance. This result is new and has a significant implication, not only in the sense discussed in the previous paragraphs regarding the theoretical epidemiological studies but also for the type of ARCH dependence found.

Indeed, although, for instance, in financial time series, the ARCH stochastic markers are long range with strongly persistent volatility, in the influenza data, those markers showed a short-range dependence that was mainly driven by the autoregressive dependence upon the squared noise. This means that the nonstationary dynamics in the variance may be linked to short-range epidemiological factors associated with disease propagation in the population, which can introduce bursts in viral transmission during an outbreak process that amplifies the dispersion and the uncertainty during those outbreak waves.

The uncovering of the ARCH effect in the influenza surveillance is new and needs further research. Since empirical studies usually focus on uncovering the chaotic dynamics and filter out the noise not studying it, which is the empirical gap that the decomposition method that we applied addresses, further studies of decomposition for influenza and other viruses are relevant to determine whether ARCH dependence is a general feature.

Although the original decomposition method that we applied here was applied to air traffic and sunspot chaos [24,25] in order to model the stochastic component, in this work, the decomposition method was adapted to a direct predictive epidemiological surveillance scheme operating on the surveillance time series data. This was achieved by testing multiple prediction algorithms using a rolling window framework that is effective when dealing with small datasets and complex nonlinear dynamics such as what we were dealing with. This brings us to the third level of contribution of this work.

There are two implications of the decomposition method for epidemiological risk analysis and surveillance. First, if we wish to extract outbreak probabilities, we need to account for both the chaotic component that affects the mean and the stochastic ARCH dynamics that affect the variance. In this way, epidemiological risk measures such as the standard deviation of the target epidemiological series and outbreak probabilities change with time, and this dynamics needs to be taken into account in epidemiological risk analysis.

The second implication for epidemiological risk analysis and surveillance regards the prediction of the time-varying mean of the process, which, in this case, results from the chaotic dynamics.

As we showed, a rolling window scheme with a machine learning unit can be used effectively to predict the mean and provide for an expectation of the target influenza surveillance series with a very high R². However, there are a few caveats, as we now explain.

For the influenza surveillance target that we worked on, which is the weekly share of positive tests for influenza in the Northern Hemisphere from January 5, 2009, to March 31, 2025, we found that, for 1-week and 2-week-ahead predictions, using the reconstructed attractor as a predictor, the adaptive local linear ridge regressor achieved R² scores of 92.11% and 85.95%, respectively. This means that the adaptive local linear approximation works well for 1 or 2 weeks; however, for 3 up to 5 weeks ahead, the ridge regressor is no longer the best performer.

In the compared models, we found the gradient boosting regressor with 100 estimators performed better than the rest for more than 2 weeks ahead, with R² scores of 81.75%, 77.59%, and 73.35% for 3-week, 4-week, and 5-week-ahead predictions, respectively.

Another relevant point is that the topological structure of the attractor contains sufficient information in the k-nearest neighbor structure to allow for up to 5-week-ahead predictions with R² scores of 78.33% (1 week), 76.16% (2 weeks), 73.88% (3 weeks), 71.51% (4 weeks), and 69.05% (5 weeks), which means that, overall, the topological structure of the reconstructed attractor supports epidemiological prediction for influenza. This has value for epidemiological planning that uses epidemiological surveillance variables such as the share of positive tests, with a high R² for up to 1 month ahead.

This study focused on a single surveillance series and region as a first exploratory step. Further studies are needed in multiple surveillance regions and different epidemiological series on influenza from a country level to a regional level. This will help achieve a deeper understanding and characterization of the influenza epidemiological dynamics regarding both the presence of chaotic dynamics and turbulent ARCH effect found in this series.

In particular, if both low-dimensional black chaos and similar ARCH dynamics are found in both hemispheres, transmission zones, and country-level data, then there is robust evidence for modeling influenza using stochastic chaos with ARCH dynamics incorporated in the noise variables.

Although black chaos is expected to be present due to the long-range intermittent process associated with viral outbreaks, the critical issue is the presence or absence of the ARCH dynamics in different surveillance series.

Similarly, for other viruses of interest, such as SARS-CoV-2, the application of the decomposition method is needed in order to uncover the type of deterministic and stochastic dynamics as well as the factors that may be present in different epidemiological data and epidemiological surveillance target variables. This is also needed to determine the implications for epidemiological surveillance, prediction, and risk analysis, since, although black chaos was also found for SARS-CoV-2, ARCH needs to be researched.