The manuscript under review presents an empirical methodology for studying stochastic chaos in epidemiological data by combining topological data analysis, topological machine learning, and nonlinear time series analysis to decompose influenza dynamics into deterministic chaotic and stochastic components down to the noise of IID random variables.

Response: Thank you for the valuable comments. The article has been revised on key points, considering the three points raised by the reviewer as explained below.

I recommend that the authors address the following shortcomings to improve the manuscript before publication.

1. The estimated largest Lyapunov exponent is approximately 0.001, which the authors characterize as “weak chaos.” However, such low values may be statistically indistinguishable from colored noise or stochastic processes with long-range dependencies. The authors should conduct additional tests using surrogate data to reliably confirm the deterministic chaotic nature of the reconstructed attractor.

Response: Surrogate tests are now applied to the original attractor and the noise-filtered attractor using both the k-nearest neighbors predictions as a noise filter, as explained above, and using a wavelet filter for comparison. Surrogate data are produced with Monte Carlo simulation, calculating bootstrap P values for both a full power law stochastic process with the same decay and for a stochastic iterative amplitude-adjusted Fourier transform–generated signal, which is more adapted than the original iterative amplitude-adjusted Fourier transform when dealing with long memory stochastic processes with different spectral features and better captures the spectrum of this series. In all of the tests, the largest Lyapunov exponents are found to be strongly statistically significant, with the noise-filtered attractors being statistically significant at the 1% level.

2. The authors should elaborate on how identifying stochastic chaos improves decision-relevant forecasting compared to established time series models in disease outbreak scenarios, including a discussion of lead times, calibration, and uncertainty quantification.

Response: Established time series forecasting models are now used in the decomposition stage, the last stage of the method, involving different families of forecasters that not only provide insights into the epidemiological dynamics but also allow us to find how different types of forecasters work for different horizons. The following models are now tested on forecasting up to 5 steps ahead (5 weeks ahead).

The linear ridge regression with penalty 0, when using the phase space reconstruction, is directly equivalent to a rolling window autoregressive process due to the fact that the delay embedding involves past values of the main series. The only issue is that since the attractor is strongly nonlinear and turbulent in the mean process and has nonstationary variance and heteroskedasticity with heavy tails even when removing the ARCH component, the autoregressive parameters are dynamical rather than fixed, and Gaussian likelihood methods are statistically incorrect for the current case.

On one step ahead (1 week) and two steps ahead (2 weeks) forecasting, the adaptive linear forecaster using linear ridge regression with penalty 0.5 is the best performer, with R² scores of 92.11% (1 week ahead) and 85.95% (2 weeks ahead) forecasting. After that point, the gradient boosting regressor model performs better.

The rolling window ridge regressor is only a local adaptive linear approximation of the nonlinear attractor, which makes its coefficients actually dynamical. This reinforces the nonlinearity of the dynamics, as explained in the article. In this case, the major dynamical coefficients, rather than being stable, fluctuate significantly, which shows the nonlinearity, but they do so in an epidemiologically meaningful way that leads to insights into the overall dynamics.

Only a rolling window is applied since a full sample estimation would defeat the purpose of having a locally adaptive test of the ability of the attractor to support adaptive forecasting for epidemiological surveillance, which is the modification done.

However, the exponentially weighted moving average (EWMA) method that is applied to capture the dynamical variance in the stochastic process is tested against the major full ARCH time series models, which allows us to test the main time series models used when dealing with nonstationary variance. It turns out that none of these models is able to capture the full process down to IID, while the rolling window ridge regressor plus EWMA model is.

On another note, since all models are using a rolling window framework with out-of-sample prediction, that is, they are trained in a rolling window and used to predict new data adapting to the attractor’s epochs and to structural shocks like the impact of COVID-19’s period on influenza, all of the reported R² scores are now forecasting scores on out-of-sample data, which did not take place in the original article in which the final R² score was a goodness-of-fit score.

The link between major epidemiological risk metrics and underlying epidemiological dynamics is addressed in greater detail in the article, namely, the SD, which in this case depends specifically on the ARCH dynamics, and the outbreak probability, which depends on both the ARCH dynamics and the chaotic process in the mean. The implications of the results are addressed in the Discussion section in relation to both the epidemiological literature on chaos in epidemiological modeling and in relation to other empirical works developed around COVID-19 and influenza. This literature is now also reviewed in detail in the Introduction section for further context.

Two appendices were added to the work, including one explaining the phase space predictive topological reconstruction using k-nearest neighbors, which comes from the original method that was initially applied to air traffic and sunspot data, and that is now reviewed in the Introduction section for more detailed context.

The role of the k-nearest neighbors algorithm, as is now explained in the first appendix, is linked to the topological structure of nonlinear dynamics and chaos, where recurrence patterns leave a predictive structure that can be captured upon appropriate embedding even without knowing the true chaotic model. These recurrence patterns are where the k-nearest neighbors algorithm is being used as a phase space reconstruction parameter search algorithm, which is more resilient to noise and different types of chaotic dynamics (including color chaos) with stochastic components. This role is now described in detail in the first appendix.

In the second appendix, we present the equations for the major ARCH models that we tested using Python’s ARCH library. In the case of the fractionally integrated generalized ARCH model, the d parameter is automatically estimated by the algorithm used.

3. The study reports high explanatory power for the full model but does not compare its performance to benchmark forecasting methods. Including a comparative analysis with traditional epidemiological or statistical models would better contextualize the added value of the topological machine learning approach and strengthen claims of methodological superiority.

Response: The method applied in this article is a 3-stage method. A significant change was applied in the last stage as a consequence of the reviewer’s comment to integrate forecasting into the method.

Indeed, in the original article, available as a preprint, we were not proposing a new method for forecasting; we were applying a recent pre-existing 3-stage method for decomposing a noisy complex attractor down to IID noise, where topological machine learning is used as a search algorithm for finding embedding parameters where there is strongest topological exploitable information linking the last phase space point to the next value of the target series.

That is, we applied a rolling window topological machine learning algorithm to find for which embedding parameters there is the highest topological structure linking a reconstructed phase space dynamic to a target time series.

This method was originally proposed and applied to air traffic data and the impact of COVID-19 on air traffic and presented at an international artificial intelligence conference; after that, it was also applied to sunspot data. A full review of this method and its previous applications has now been introduced in the revised article for better clarity.

The original method has the following three stages:

The reported final R² that we presented for the preprint version was a goodness-of-fit metric rather than an out-of-sample forecasting score because the residuals’ stochastic process was researched and estimated using the full sample.

Now, while this is what we applied in the original preprint article, we made a few changes in the revised article, reflecting the reviewer’s comment, which opened up the possibility of adapting the method for epidemiological forecasting, thus adding value for epidemiological surveillance and risk analysis.

This is a different approach, implying the need for a change at the decomposition stage (third stage of the method). In this instance, rather than using the filtered residuals with the k-nearest neighbors predictions used as an adaptive topological forward-looking noise filter and fitting a time series model to the resulting noise down to IID components, we changed the third stage of the model to test different forecasting schemes that can be implemented using the reconstructed attractor and the resulting residuals to capture the full process down to IID noise in a way that allows for forecasting and supports risk analysis and epidemiological surveillance rather than just epidemiological characterization, which was the original intention of the article.

In this way, while the first two stages of the method are implemented as in the original method, the third stage has been changed to incorporate a forecasting-based decomposition that provides value for epidemiological surveillance since the decomposition stage directly targets surveillance and risk analysis while at the same time covering the modeling component.

To perform a decomposition that can be turned into an epidemiological surveillance scheme and to investigate the degree to which the reconstructed attractor can support epidemiological surveillance without losing the decomposition and modeling dimension, that is, the point of the last stage was the critical change that led to a more extensive revision of the empirical section at the last stage of the method.

The revised article in the Results section now features four subsections: the first is the characterization of the surveillance series with the signal spectral analysis, and subsections three to four correspond to the three stages of the method.

It is at the last stage that we implemented the switch from the original method, which was developed for “dynamics characterization,” to a “dynamics characterization and forecasting” approach.

In this case, we test different prediction algorithms on the reconstructed attractor using the same rolling window as the one used for the adaptive k-nearest neighbors model and use the k-nearest neighbors as a benchmark to get a baseline on topological predictability on multiple forecasting horizons.

The main point at the first part of the third stage, before decomposition, is now to find the degree to which the attractor supports prediction on different horizons and which models work best for which horizons, covering different major families of models.

These models are, thus, first tested on prediction on a 1 week ahead prediction horizon using the reconstructed attractor, and then they are also tested on multiple weeks ahead prediction to evaluate how the attractor can support prediction on different horizons and which algorithms should be used.

It becomes evident from the results that there is no single best algorithm that works for all horizons. This is due to the nature of the attractor; thus while for short horizons we get best results from one algorithm, for longer horizons, we have to switch to another algorithm for forecasting and epidemiological surveillance.

After addressing this part of the forecasting issue, which is still within the context of characterizing the nature of the chaotic dynamics but at the same time has direct implications for epidemiological surveillance, we proceed to the decomposition process using the best performer in the one period ahead forecasting and using the residuals to address the stochastic component.

In this last stage, again, we introduce a change to the original method and article to best fit the forecasting problem to get a better fit into epidemiological surveillance rather than just epidemiological modeling. Therefore, we no longer estimate a model on the full sample residuals, which is what is done in the original method that is aimed at supporting stochastic model testing on a filtered noise term rather than forecasting.

Instead, we now evaluate a dynamical method that does not require the full sample and that can be integrated into the adaptive framework of the article and of the resulting epidemiological surveillance process.

Specifically, we still find nonstationary variance with ARCH signatures in the residuals of the best performer in the one-period forecasting, so we implement and test an EWMA ARCH model, which is better suited for the forecasting system that is used in the rolling window machine learning framework. Additionally, we test whether the EWMA algorithm is capable of leading to IID noise and of fully capturing the nonstationary variance component consistent with the forecasting scheme, which it does.

We also compare the EWMA’s performance with major ARCH models, showing that neither of these models leads to IID noise. With this modification, the method used not only allows us to decompose the attractor down to IID noise (original purpose of the method) but also allows us to draw direct insights into epidemiological surveillance and risk analysis in a more effective way than in the original preprint version.

The last stage of the model is thus renamed from decomposition to predictive decomposition to highlight the change to the original method. The references to the original three-stage method are presented, and this change to the method is now addressed in the revised article.

It is important to stress that the comparison of the prediction algorithms tested in the third stage of the method cannot be moved to the first stage of the method with the algorithms used for deciding the phase space embedding. Since the first stage of the method is not about a comparison of prediction algorithms, instead it is about finding which embedding captures the strongest topological structure linking the phase space trajectory to the time series, so a topological machine learning unit has to be used in the first stage of the method. The employment of k-nearest neighbors for this is effective since it is a noise-resilient algorithm, making it more appropriate when dealing with stochastic chaos. This explanation is now integrated for greater clarity into the Methods section as well as in Appendix A.

In the last stage of the method, different algorithms can indeed be tested, since this is a decomposition stage rather than an embedding selection stage. So the inclusion of these different algorithms makes the predictive decomposition and, in this way, adds value from an epidemiological surveillance standpoint, which now constitutes the third stage of the modified method.