A Koopman operator-based prediction algorithm and its application to COVID-19 pandemic and influenza cases … – Nature.com

We apply our algorithms to a few case studies in epidemiology: Influence epidemics and COVID-19. We do emphasize that the techniques are general and can be applied to any system that experience a drastic change in its fundamental behavior.

As first example for showing our prediction methodology, we use the set of data associated with influenza epidemics. Clearly, not driven by an underlying deterministic dynamical system, the influenza time series exhibits substantial regularity in that it occurs typically during the winter months, thus enabling coarse-grained prediction of the type we will see a very small number of cases of influenza occurring in summer months. However, predicting the number of influenza cases accurately is a notoriously hard problem20, exacerbated by the possibility that a vaccine designed in a particular year does not effectively protect against infection. Moreover, the H1N1 pandemic that occurred in 2009 is an example of a Black Swan event.

The World Health Organizations FluNet is a global web-based tool for influenza virological surveillance. FluNet makes publicly available data on the number of specimens with the detected influenza viruses of type A and type B. The data have been collected from different countries, starting with the year 1997, and are updated weekly by the National Influenza Centers (NICs) of the Global Influenza Surveillance and Response System (GISRS) and other national influenza reference laboratories, collaborating actively with GISRS. We use the weekly reported data for different countries, which consist of the number of received specimens in the laboratories, the distribution of the number of specimens with confirmed viruses of type A.

The Koopman Mode Decomposition was used in the context of analyzing the dynamics of the flu epidemic from differentGoogle Fludata in21. We remark that the authors of that paper have not attempted prediction, and have analyzed only stationary modese.g. the yearly cycles, thus making the papers goals quite different from the nonstationary prediction pursued here.

We first compare the global and the local prediction algorithms. The KMD is computed using active windows of size ({textsf{w}}= 312), and the (208 times 104) HankelTakens matrices. In Fig.reff1a, we show the performances of both algorithms, using the learning data from the window April 2003April 2009 (shadowed rectangle). In the global prediction algorithm the dynamics is predicted for 104 weeks ahead. The first type of failure in the global prediction algorithm and forecasting appears after the Black Swan event occurred in the years 2009 and 2010. This is recognized by the algorithm, so that it adapts by using the smallest learning span and, with this strategy, it allows for reasonably accurate forecasting, at least for shorter lead times. This data, in addition to those from Supplementary Information section S2.4 show the benefits of monitoring the prediction error and switching to local prediction. The initial HankelTakens matrix is (3times 2), and the threshold for the local prediction relative error in Supplementary Information Algorithm S4 is 0.005.

Influenza data (USA). (a) The data are collected in the window April 2003April 2009 (shadowed rectangle) and then the dynamics is predicted for 104 weeks ahead. The local prediction algorithm recovers the prediction capability by forgetting the old data and using narrower learning windows. The local prediction algorithm delivers prediction for one week ahead. (b) The active window (shadowed rectangle) is July 2004July 2010, and the dynamics is predicted for 104 weeks ahead. The global prediction fails due to the Black Swan data in the learning window. (Some predicted values were even negative; those were replaced with zeros.) The global prediction algorithm recovers after the retouching the Black Swan event data, which allows for using big learning window. Compare with positions of the corresponding colored rectangles in Fig.2.

Next, we introduce an approach that robustifies the global algorithm in the presence of disturbances in the data, including the missing data scenario. We use the data window July 2004July 2010, which contains a Black Swan event in the period 20092010. As shown in Fig.1b, the learned KMD failed to predict the future following the active training window. This is expected because the perturbation caused by the Black Swan event resulted in the computed Ritz pairs that deviated from the precedent ones (from a learning window before disturbance), and, moreover, with most of them having large residuals. This can be seen as a second type of failure in the global prediction.

The proposed Black Swan event detecting device, built in the prediction algorithm (see Supplementary Information Algorithm S3), checks for this anomalous behaviour of the Ritz values and pinpoints the problematic subinterval. Then, the algorithm replaces the corresponding supplied data with the values obtained as predictions based on the time interval preceding the Black Swan event. Figure1b shows that such a retouching of the disturbance allows for a reasonable global prediction.

Note that in a realistic situation, global predictions of this kind will trigger response from authorities and therefore prevent its own accuracy and induce loss of confidence, whereas local prediction mechanisms need to be deployed again.

The real and imaginary parts of Ritz values with residuals bellow (eta _r=0.075) for sliding active windows. The color intensity of eigenvalues indicates the amplitudes of the corresponding modes. Pink rectangles mark ends of training windows with no acceptable Ritz values. Note how the unstable eigenvalues ((Re (lambda )>0)) impact the prediction performance, and how the retouching moves them towards neutral/stablethis is shown in the yellow rectangle in panels (a) and (c). Also influenced by the disturbance are the eigenvalues in the light blue rectangles in panels (a), (b); retouching moves the real parts of eigenvalues towards neutral/stable and rearranges them in a lattice-like structure22, as shown in panels (c), (d). Compare with Fig.1b.

We now discuss the effect of the Black Swan event and its retouching to the computed eigenvalues and eigenvectors. We have observed that, as soon as a disturbance starts entering the training windows, the Ritz values start exhibiting atypical behavior, e.g. moving deeper into the right half plane (i.e. becoming more unstable), and having larger residuals because the training data no longer represent the Krylov sequence of the underlying Koopman operator.

This is illustrated in the panels (a) and (b) in Fig.2, which show, for the sliding training windows, the real and the imaginary parts of those eigenvalues for which the residuals of the associated eigenvectors are smaller than (eta _r=0.075). Note the absence of such eigenvalues in time intervals that contain the disturbance caused by the Black Swan event.

On the other hand, the retouching technique that repairs the distorted training data restores the intrinsic dynamics over the entire training window. The distribution of the relevant eigenvalues becomes more consistent, and the prediction error decreases, see panels (c) and (d) in Fig.2, and in Supplementary Information Figure S16.

Our proposed retouching procedure relies on detecting anomalous behavior of the Ritz values; a simple strategy of monitoring the spectral radius of active windows (absolutely largest Ritz value extracted from the data in that window) is outlined in Supplementary Information. Note that this can also be used as a litmus test for switching to the local prediction algorithm. In Supplementary Information, we provide further examples, with the influenza data, that confirm the usefulness of the retouching procedure. In general, this procedure can also be adapted to the situation when the algorithm receives a signal that the incoming data is missing or corrupted.

The second set of data we consider is that associated with the ongoing COVID-19 pandemic. Because the virus is new, the whole event is, in a sense, a Black Swan. However, as we show below, the prediction approach advanced here is capable of adjusting quickly to the new incoming, potentially sparse data and is robust to inaccurate reporting of cases.

At the beginning of the spread of COVID-19, we have witnessed at moments rather chaotic situation in gaining the knowledge on the new virus and the disease. The development of COVID-19 diagnostic tests made tracking and modeling feasible, but with many caveats: the data itself is clearly not ideal, as it depends on the reliability of the tests, testing policies in different countries (triage, number of tests, reporting intervals, reduced testing during the weekends), contact tracing strategies, using surveillance technology, credit card usage and phone contacts tracking, the number of asymptomatic transmissions etc. Many different and unpredictable exogenous factors can distort it. So, for instance the authors of23 comment at https://ourworldindata.org/coronavirus-testing that e.g. The Netherlands, for instance, makes it clear that not all labs were included in national estimates from the start. As new labs get included, their past cumulative total gets added to the day they begin reporting, creating spikes in the time series. For a prediction algorithm, this creates a Black Swan event that may severely impair prediction skills, see sectionRetouching the Black Swan event data.

This poses challenging problems to the compartmental type models of (SIR, SEIR) which in order to be useful in practice have to be coupled with data assimilation to keep adjusting the key parameters, see e.g.24. Our technique of retouching (sectionRetouching the Black Swan event data) can in fact be used to assist data assimilation by detecting Black Swan disturbance and thus to avoid assimilating disturbance as normal.

In the KMD based framework, the changes in the dynamics are automatically assimilated on-the-fly by recomputing the KMD using new (larger or shifted) data snapshot windows. This is different from the compartmental type models of infectious diseases, most notably in the fact that the procedure presented here does not assume any model and, moreover, that it is entirely oblivious to the nature of the underlying process.

As a first numerical example, we use the reported cumulative daily cases in European countries. In Supplementary Information section S1.5, we use this data for a detailed worked example that shows all technical details of the method. This is a good test case for the methodusing the data from different countries in the same vector observable poses an additional difficulty for a data driven revealing of the dynamics, because the countries independently and in an uncoordinated manner impose different restrictions, thus changing the dynamics on local levels. For instance, at the time of writing these lines, a new and seemingly more infectious strain of the virus circulating in some parts of London and in south of England prompted the UK government to impose full lockdown measures in some parts of the United Kingdom. Many European countries reacted sharply and immediately suspended the air traffic with the UK.

In the first numerical experiment, we use two datasets from the time period February 29 to November 19. and consider separately two sets of countries: Germany, France and the UK in the first, and Germany, France, UK, Denmark, Slovenia, Czechia, Slovakia and Austria in the second. The results for a particular prediction interval are given in Figs.3 and 4. For more examples and discussion how the prediction accuracy depends on the Government Response Stringency Index (GRSI25,26) see Supplementary Information section S1.5.

Prediction of COVID-19 cases (35 days ahead, starting July 11) for Germany, France and United Kingdom. Left panel: The HankelTakens matrix ({mathbb {H}}) is (282 times 172), the learning data consists of ({textbf{h}}_{1:40}). The KMD uses 39 modes. Middle panel: The matrix ({mathbb {H}}) is (363 times 145), the learning data is ({textbf{h}}_{1:13}). The KMD uses 12 modes. Right panel: The KoopmanRitz values corresponding to the first (magenta circles) and the middle (blue plusses) panel. Note how the three rightmost values nearly match.

Prediction errors and KMD spectrum of COVID-19 cases (28 days ahead, starting July 11) for Germany, France, United Kingdom, Denmark, Slovenia, Czechia, Slovakia and Austria. Left panel: The HankelTakens matrix ({mathbb {H}}) is (752 times 172), the learning data consists of ({textbf{h}}_{1:40}). The KMD uses 39 modes. Middle panel: The matrix ({mathbb {H}}) is (968 times 145), the learning data is ({textbf{h}}_{1:13}). The KMD uses 12 modes. Right panel: The KoopmanRitz values corresponding to the first two computations in Fig.3 (magenta circles and blue pluses, respectively) and the the first two panels in this Figure (orange x-es and cyan squares, respectively). Note how the corresponding KoopmanRitz values nearly match for all cases considered.

In the above examples, the number of the computed modes was equal to the dimension of the subspace of spanned by the training snapshots, so that the KMD of the snapshots themselves was accurate up to the errors of the finite precision arithmetic. In general, that will not be the case, and the computed modes will span only a portion the training subspace, meaning that the KMD of the snapshots might have larger representation error. (Here we refer the reader to Supplementary Information section S1.3, where all technical details are given.) This fact has a negative impact to the extrapolation forward in time and the problem can be mitigated by giving more importance to reconstruction of more recent weights. This is illustrated in Figs.5 and6, where the observables are the raw data (reported cases) for Germany, extended by a two additional sequence of filtered (smoothened) values.

Prediction experiment with data from Germany. Left panel: the computed residuals for the computed 102 Koopman Ritz pairs (extracted from a subspace spanned by 132 snapshots ({textbf{h}}_{1:132})). Note that all residuals are small. The corresponding Ritz values are shown in the first panel in Fig.6. Middle panel: KMD reconstruction error for ({textbf{h}}_{1:132}) and the error in the predicted values ({textbf{h}}_{133:160}) (encircled with ({circ })). The reconstruction is based on the coefficients ((alpha _j)_{j=1}^r=mathrm {argmin }_{alpha _j}sum _{k} Vert {textbf{h}}_k - sum _{j=1}^{r} lambda _j^{k}alpha _j {textbf{v}}_jVert _2^2). Right panel: Prediction errors for the period October 11November 7.

Prediction experiment with DS3 with data from Germany. Left panel: the computed 102 Koopman Ritz values (extracted from a subspace spanned by 132 snapshots ({textbf{h}}_{1:132})). The corresponding residuals are shown in the first panel in Fig.5. Middle panel: KMD reconstruction error for ({textbf{h}}_{1:132}) and the error in the predicted values ({textbf{h}}_{133:160}) (encircled with ({circ })). The reconstruction is based on the coefficients ((alpha _j)_{j=1}^r=mathrm {argmin }_{alpha _j}sum _{k} w_k^2 Vert {textbf{h}}_k - sum _{j=1}^{r} lambda _j^{k}alpha _j {textbf{v}}_jVert _2^2). Right panel: Prediction errors for the period October 11 November 7. Compare with the third graph in Fig.5.

The figures illustrate an important point in prediction methodology, that we emphasized in the introduction: a longer dataset and a better data reconstruction ability (i.e. interpolation) does not necessarily lead to better prediction. Namely, weighting more recent data more heavily produces better prediction results. This was already observed in27 for the case of traffic dynamics, and the method we present here can be used to optimize the prediction ability.

We have deployed the algorithm to assess the global and United States evolution of the COVID-19 pandemic. The evolution of the virus is rapid, and Black Swans in the sense of new cases in regions not previously affected appear with high frequency. Despite that, the Koopman Mode Decomposition based algorithm performed well.

In Fig.7a we show the worldwide forecast number of confirmed cases produced by the algorithm for November 13th, 2020. The forecasts were generated by utilizing the previous three days of data to forecast the next three days of data for regions with higher than 100 cases reported. The bubbles in Fig.7a are color coded according to their relative percent error. As can be observed, a majority of the forecasts fell below 15% error. The highest relative error for November 13th, 2020 was 19.8% which resulted from an absolute error of 196 cases. The mean relative percent error, produced by averaging across all locations, is 1.8% with a standard deviation of 3.36% for November 13th, 2020. Overall, the number of confirmed cases are predicted accurately and since the forecasts were available between one to three days ahead of time, local authorities could very well utilize our forecasts to focus testing and prevention measures in hot-spot areas that will experience the highest growth.

A video demonstrating the worldwide forecasts for March 25, 2020November 29, 2020 is provided in the Supplementary Information online (Fig.7a is a snapshot from that video). Lastly, it is well known that the ability to test people for the virus increased throughout the development of the pandemic and thus resulted in changes in the dynamics of reported cases. Although it is impossible for a data-driven algorithm to account for changes due to external factors, such as increased testing capabilities, it is important that the algorithm be able to adjust and relearn the new dynamics. For this reason, we encourage the reader to reference the video and note that although periods of inaccuracy due to black swan events occur, the algorithm is always able to stabilize and recover. In contrast, since this is at times a rapidly (exponentially) growing set of data, methods like naive persistence forecast do poorly.

In Fig.7b, c we show the performance of the prediction for the cumulative data for the US in March-April 2020. It is of interest to note that the global curve is obtained as a sum of local predictions shown in Fig.7a, rather than as a separate algorithm on the global data. Again, the performance of the algorithm on this nonstationary data is good.

Prediction of confirmed COVID-19 cases utilizing the publicly available COVID-19 data repository provided by Johns Hopkins. The true data ranges between March 22nd, 2020 and November 29th, 2020. We utilize the last three days of data to forecast the following three days of data. (a) Predicted conditions and prediction error worldwide on November 13. The widths of the bubbles represent the number of cases in a region; only regions with more that 100 cases are used and the bubbles are colored according to their relative percent error. (b) Comparison of true and forecast data for cumulative confirmed cases in the US for April to December 2020. The cumulative forecasts shown here were obtained by summing the forecasts of the individual locations, indicating that the region specific forecasts were sufficiently accurate for tracking the cumulative dynamics of the virus in the US. (c) Percent error for the forecasts of the cumulative confirmed cases in the US. On average the percent error is less than 5 percent and although spikes occur, which could be due to changes in testing availability, the algorithm adjusts and the error stabilizes within a short amount of time. Furthermore, Johns Hopkins provided data for around 1787 locations around the United States and we produced forecasts for each of those locations.

Read more here:

A Koopman operator-based prediction algorithm and its application to COVID-19 pandemic and influenza cases ... - Nature.com