^{1}

^{2}

^{*}

^{1}

^{2}

^{1}

^{2}

Edited by: Ulrich Parlitz, Max-Planck-Institut für Dynamik und Selbstorganisation, Germany

Reviewed by: Hiromichi Suetani, Oita University, Japan; Xin Tong, National University of Singapore, Singapore

This article was submitted to Dynamical Systems, a section of the journal Frontiers in Applied Mathematics and Statistics

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Data assimilation permits to compute optimal forecasts in high-dimensional systems as, e.g., in weather forecasting. Typically such forecasts are spatially distributed time series of system variables. We hypothesize that such forecasts are not optimal if the major interest does not lie in the temporal evolution of system variables but in time series composites or features. For instance, in neuroscience spectral features of neural activity are the primary functional elements. The present work proposes a data assimilation framework for forecasts of time-frequency distributions. The framework comprises the ensemble Kalman filter and a detailed statistical ensemble verification. The performance of the framework is evaluated for a simulated FitzHugh-Nagumo model, various measurement noise levels and for

Understanding the dynamics of natural complex systems is one of the great challenges in science. Various research domains have developed optimized analytical methods, computational techniques or conceptual frameworks to gain deeper insight into the underlying mechanisms of complex systems. In the last decades, more and more interdisciplinary research attracted attention building bridges between research domains by applying methodologies outside of domains. These cross-disciplinary techniques fertilize research domains and shed new light on their underlying properties. A prominent example is the mathematical domain of dynamical systems theory that traditionally is applied in physics and engineering and that has been applied very successfully in biology and neuroscience. For instance, taking a closer look at the spatio-temporal nonlinear dynamics of neural populations has allowed to identify epilepsy as a so-called dynamical disease [

Weather forecasts are an everyday service provided by national and regional weather services that allows to plan business processes as well as private activity and serves as a warning system for extreme weather situations, such as floods or thunderstorms. Weather forecast is also an important research domain in meteorology that has been developed successfully in the last decades improving the forecasts for both global phenomena and local weather situations. In detail, todays weather services employ highly tuned and optimized meteorological models and data processing techniques to compute reliable forecasts. Specifically the combination of an efficient model and measured meteorological data enables researchers to provide various types of predictions, such as the probability of rain or the expected temperature in certain local regions. This optimal combination of model and data is achieved by data assimilation [

In other research domains, prediction methods are rare but highly requested. For instance, the prediction of epileptic seizures [

Most recent data assimilation studies apply the unscented Kalman filter [

The work is structured as follows. The Methods section introduces the model, simulated observations, the ensemble Kalman filter and the verification metrics applied. The subsequent section shows obtained results for

Single biological neurons may exhibit various types of activity, such as no spike discharge, discharge of single spikes, regular spike discharge, or spike burst discharges. These activity modes can be described by high-dimensional dynamical models. A more simple model is the FitzHugh-Nagumo model [

with membrane potential

with maximum time _{n} and _{n} over time results in a shift of oscillation frequency of the system, i.e., from larger to smaller frequencies. Such a non-homogeneous temporal rhythm is well-known in neuroscience, e.g., in the presence of anesthetic drugs [

leading to a single oscillation frequency. We point out that τ_{n}(_{f} and _{n}(_{f} and both models converge to each other for

The model integration over time uses a time step of 0.01 and every 50 steps a sample is written out running the integration over 5·10^{4} steps in total. Initial conditions are ^{t}. After numerical integration, we re-scaled the unit-less time by α_{max} = 1s. This sets the number of data points to

To reveal non-stationary cyclic dynamics, we analyze the time-frequency distribution of data with spectral density _{k}, ν_{m}),

applied uses a mother wavelet Ψ with central frequency _{c} = 8 and the time-frequency distribution has a frequency resolution of Δν = 0.5 Hz in the range ν ∈ [5 Hz;20 Hz]. The parameter _{c}/ν is the scale that depends on the pseudo-frequency ν. By the choice of the central frequency _{c}, the mother wavelet has a width of 4 periods of the respective frequency. This aspect is important to re-call while interpreting temporal borders of time-frequency distributions. For instance, at a frequency of 15 Hz border disturbances occur in a window of 0.26 s from the initial and final time instant.

Figure

Noise-free dynamics of the Fitzhugh-Nagumo model and corresponding observations.

To relate model variables to observations, data assimilation introduces the notion of a measurement operator

The system dynamics can be observed in various ways and the observation operator is chosen correspondingly. Measurements directly in the system are called

The present study considers scalar

where ξ(

Noisy

From Equation (4), one reads off the observation operator

with ^{t} ∈ ℜ^{2}.

For comparison, we also consider nonlocal observations with the observation operator

yielding

for the same noise levels κ as above. Figure

Nonlocal and speed observations at various measurement noise levels κ.

As already stated, the aim of the present work is to introduce the idea to forecast temporal features. As a further step in this direction, let us consider temporal changes of the signal evolution, i.e., the speed of the system. To this end the definition of the observation operator

with

yielding

for two noise levels κ = 0.0 and κ = 0.02. Numerically, the derivative _{n})−_{n−1}) at time instance _{n}. Figure

One of the major aims of data assimilation techniques is the optimal fit of model dynamics to observed data. Here, we introduce the major idea with a focus on the 2-dimensional model (1) and the scalar observation. Observations

To merge observation _{b}(_{a} minimizes the cost function

i.e., the solution is the minimum of the cost function _{a} is called the analysis, ^{2}. However, typically, one does not know the true observation error κ and

This is the major result of the 3DVar technique for scalar observations [

Conversely, if the covariance error matrix

with the ensemble mean

and in observation space

with _{a, b} = _{a, b}. Since ^{t} and

stating that the analysis equivalent in observation space _{a} is always closer to the observation as the background equivalent in observation space _{b}.

The ensemble transform Kalman filter (ETKF) [

with the ensemble space coordinates ^{L}. Re-considering the optimization scheme (7) in this space

with ^{L×L}. Then the analysis ensemble mean _{a} reads

The analysis ensemble members can be calculated by

with

corresponding to (10) and with ^{l} ∈ ℜ^{L}. Defining the matrix ^{L×L} with columns ^{l}, the ansatz ^{t}, and Equation (15) yields ^{t}, the orthogonal matrix

where

Equation (7) implies that all states, observations, covariances and operators are instantaneous. Extensions of this formulation are known, e.g., as the 4D-ENKF or the 4DVar [

In each analysis step, the analysis equivalent in observation space _{a} moves away from the model background state _{b} closer to the observation ^{t}. However, the model has errors which are not completely reflected by the state estimate error covariance matrix ^{t}, since this is calculated based on an ensemble of model forecasts with the same simulated model equations. To take care of the model error and draw the analysis closer to the background state, typically one enhances the ensemble spread by inflation.

For

Putting together models and data assimilation, the model evolution is controlled by observed data optimizing the initial state of the model iteration. Our data assimilation cycle starts with initial conditions from which the model evolves during the sampling interval. The model state after one sampling interval Δ_{b}. The subsequent data assimilation step estimates the analysis state _{a} that represents the initial state for the next model evolution step. In other words, data assimilation tunes the initial state for the model evolution after each sampling interval. Using the ETKF, this cycling is applied for all ensemble members which obey the model evolution and whose analysis state is computed in each data assimilation step. Initial ensemble member model states were _{1}, η_{2} in the range η_{1}, η_{2} ∈ [0;1].

The aim of the present work is to show how optimal forecasting can be done. Free ensemble forecasts are model evolutions over a time typically longer than the sampling time. This forecast time is called lead time. The initial state of the free forecasts are the analysis model states determined by data assimilation.

In the present work, we are interested in forecasts at every sample time instant. To this end we compute the model activity at a certain lead time. This forecast is computed for all ensemble members what renders it an ensemble prediction. The forecasts are solutions of the model _{a} with initial analysis state _{a} at time _{a} and lead time _{a}. To compare them to observations, forecasts are mapped to observation space yielding model equivalents

Later sections show free forecasts ^{f}(

Naturally, one expects that the forecasts diverge from observations with longer lead times but the question is which forecasts can still be trusted, i.e., are realistic. Essentially we ask the question how one can verify the forecasts. To this end, various metrics and scores have been developed [

To estimate the deviation of forecast ensemble members _{n},

the root-mean square error

and the ensemble spread

For scalar observations and corresponding forecasts, i.e., temporal time series,

The time-frequency distribution represents the spectral power distribution _{1}(_{k}, ν), _{2}(_{k}, ν) at time instance _{k}. A corresponding well-known distance measures is the time-averaged Itakura-Saito distance (ISD) [

This distance measure is not symmetric in the spectral distributions and hence not a metric. As an alternative, one may also consider the log-spectral distance (LSD) [

which has the advantage that it is symmetric in the distributions. In both latter measures _{obs} and _{fc} are the power spectra of observations and forecasts, respectively.

As pointed out above, we hypothesize that spectral features extracted from forecasts can be predicted in a better or more precise way than forecasts themselves. Since measurement noise plays an important role in experimental data, we evaluate predictions for medium and large noise levels κ compared to κ = 0. The

reflects the deviation of forecast errors at medium and large noise levels from noiseless forecasts. For SS = 0, forecasts have identical rmse and SS < 0 (SS > 0 ) reflects larger (smaller) rmse, i.e., worse (better) forecasts. The skill score SS is less sensitive to the bias as the rmse, and that also plays an important role in the evaluation of forecasts (similarly to the standard deviation). However, for small bias SS > 0 is a strong indication of improved forecasts.

According to Equation (10), the ensemble is supposed to describe well the model error. The ensemble spread represents the variability of the model and an optimal ensemble stipulates spread = rmse [

quantifies this relation. If

A representative forecast ensemble has the same distribution as the observations. This can be quantified by computing the rank of an observation in a forecast ensemble [

with the gamma-function Γ(^{2} permits to estimate the function parameters by

The derived β−score [

equals 0 for a uniform distribution and β_{c} > 0 (β_{c} < 0) reflects the ensemble overestimation (underestimation) of the model uncertainty for an inverse U-shaped (U-shaped) distribution. In addition, the β-bias [

quantifies the skewness of the rank distribution and β_{b} = 0 reflects symmetric distributions. β−bias values β_{b} > 0 (β_{b} < 0) reflect a weight to lower (higher) ranks and the majority of ensemble members is larger (smaller) than observations.

At first, we consider

To start, we consider

To illustrate the ensemble evolution, Figure

Illustration of the temporal evolution of ensemble spread in observation space.

Now let us turn to the forecasts. In the data assimilation cycle, after one model step and hence one sampling time interval, the analysis is computed and initializes the phase space trajectory of the model evolution for the subsequent model step. In free forecasts _{a} initialized by the analysis at each time instant _{a}. Figures

The time-frequency distribution of the observations and forecast equivalents is shown in Figures

To quantify the differences between forecasts and observations detected by visual inspection in section 3.2, we compute the forecast departure statistics subject to the lead time. Figure

Ensemble verification metrices and scores with respect to lead times.

The ensemble spread decreases with lead time in both time series data and time-frequency data to values smaller than the rmse. This yields a decreasing spread-skill relation where

The reliability of the ensemble forecasts can be evaluated by rank histograms, i.e., the β−score β_{c} and β−bias β_{b}. Figure _{c} decreases from positive to negative values both for time series and time-frequency distribution data. This reveals an underestimation of the model uncertainty. The β−bias remains positive-definite for time series data whereas β_{b} of time-frequency distribution data decreases from positive to negative values. This result reveals that the majority of ensemble members are larger than the time series observations and smaller than the time-frequency spectral power observations.

To understand better why the ensemble spread shrinks at large lead time, Figure

Phase space dynamics for short and long lead times. The ensemble mean of forecasts and the true data are color-coded in red and black, respectively. The blue-coded points represent the false model data.

To understand how specific the gained results from

Time series and time-frequency distributions of free ensemble mean forecasts compared to nonlocal observations.

To understand this, we take a closer look at the forecast time series at _{0} = 12.5 Hz since then _{0} is exactly one period of this oscillation. This fixed phase relation is observed in Figure _{0} = 1/

The departure statistics between forecasts and observations resembles the findings for

The departure metrics Bias, rmse and spread of forecasts of nonlocal observations, the corresponding skill score SS and spread-skill ratio _{c} and β_{b}

These results are in good accordance to the rank histogram features β_{c} and β_{b} seen in Figure _{c} > 0 reflecting an overestimation of the spread, otherwise β_{c} < 0 reflecting a too small ensemble spread. This holds true for all data types and all noise levels. The β−bias is similar to Figure

Summarizing, the ensemble varies much with the lead time what indicates a fundamental problem in the ensemble forecast.

Spectral power takes into account data at several time instances. Since to our knowledge Kalman filters have not been developed yet for observation operators nonlocal in time, we take a first step and consider speed observations subjected to two noise levels. Figure

Observations, first guess and analysis for speed data.

To improve the assimilation cycle, we diminish the observation error to

The forecasts in Figure

Time series of observations and forecasts (model equivalent in observation space) at two lead times and two noise levels.

This can be quantified by departure statistics metrics as shown in Figure

The departure metrics bias, rmse, and spread of forecasts to speed observations and the corresponding skill score SS and spread-skill ratio

These results are in good accordance to the rank histogram features β_{c} and β_{b} seen in Figure _{c} for all lead times reflects the underestimation of the spread and the β−bias β_{b}≈0 indicates that this underestimation is present for all forecast values. This holds true for both noise levels.

Rank statistics β_{c} and β_{b} for the ensemble in the presence of speed observations. The color-coding is identical to Figure

Since the rmse is not an optimal measure to quantify the difference between time-frequency distributions, we compute more advanced measures specific for power spectra. The Itakura-Saito distance (ISD) and the log-spectral distance (LSD) increase with the lead time for

Improved verification statistics of

Time frequency distributions appear to represent instantaneous spectral power. However, the spectral power distributions at subsequent time instances are strongly correlated dependent on the frequency. The correlation length is τ = 4/

The present work applies well-established techniques known in meteorology to find out whether they can be useful to forecast spectral features in other science domains where spectral dynamics plays an important role, such as in neuroscience. For

Since time-frequency distributions show time-variant spectral power, it is necessary to verify forecasts by spectral power-specific measures and take care of spectral power-specific artifacts, cf. Figure

Conversely, speed observations consider the dynamical evolution of the system and are a very first approximation to a direct spectral feature. This is true since speed observations do not take into account the system state and observation at a single time instance only. Future work will extend this approach to a larger time window what allows to compute the power spectrum that can be mapped to a single time instance. Since generalizations or differential operators are integral operators [

Since spectral feature forecasts are sensitive to certain frequencies, they are sensitive to lead time-observation frequency resonances. Such resonances seem to improve the forecast although these resonances are artifacts. To our best knowledge, the current work is the first to uncover these resonances that may play an important role in the interpretation of forecasts.

The ensemble data assimilation cycle involves several modern techniques, such as multiplicative and additive covariance inflation that well improves the forecasts. As a disadvantage, the spread for short lead times is too large. Future work will improve the ensemble statistics by adaptive inflation factors [

The ensemble Kalman filter applied is one possible technique to gain forecasts. Other modern powerful techniques are the variational methods 3D- and 4D-Var [

Eventually, the present study considers a specific model system that exhibits a single time scale due to a single oscillation frequency. However, natural complex systems exhibit multiple time scales what may render the Kalman filter less effective and the superiority of the time-frequency data less obvious. In the future, it will be an important task to extend the present work to multi-scale Kalman filters [

AH conceived the study and performed all simulations. AH and RP planned the manuscript structure and have written the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors would like to thank Felix Fundel, Michael Denhardt, and Andreas Rhodin for valuable discussions.