Inferring the connectivity of coupled chaotic oscillators using Kalman filtering

Inferring the interactions between coupled oscillators is a significant open problem in complexity science, with multiple interdisciplinary applications. While the Kalman filter (KF) technique is a well-known tool, widely used for data assimilation and parameter estimation, to the best of our knowledge, it has not yet been used for inferring the connectivity of coupled chaotic oscillators. Here we demonstrate that KF allows reconstructing the interaction topology and the coupling strength of a network of mutually coupled Rössler-like chaotic oscillators. We show that the connectivity can be inferred by considering only the observed dynamics of a single variable of the three that define the phase space of each oscillator. We also show that both the coupling strength and the network architecture can be inferred even when the oscillators are close to synchronization. Simulation results are provided to show the effectiveness and applicability of the proposed method.

www.nature.com/scientificreports/ is different from zero, the oscillators are connected; otherwise, they are not. Here, as a first step, we show that KF can indeed infer the coupling strengths in a network of Rössler-like chaotic oscillators. To do this, we numerically simulate Rössler-like chaotic oscillators as described in Sevilla-Escoboza et al. 39 , and use the so-called Unscented Kalman Filter (UKF) 40 implemented in Matlab.
Our motivation to use the expression of the Rössler-like chaotic oscillators given in 39 is the existence of a large and freely available experimental data set (available for Ref. 39 and for Ref. 41 ). However, as we will discuss later, the temporal filtering performed by the experimental data acquisition system prevents obtaining a good network reconstruction. Therefore, here, as a proof-of-concept demonstration, we limit our study to consider synthetically-generated time series. Future work will attempt to incorporate the data acquisition process into the method presented here, to try to reconstruct the network topology using the available experimental data.
In this study, we know the details of the governing equations of each Rössler-like oscillator and leave their coupling strengths as the only unknown parameters. We first test the algorithm with an isolated oscillator and show that, from the analysis of one of the variables, it is possible to retrieve the evolution of the other two variables which describe the system. Next, we consider two mutually coupled oscillators. For these two elements, first, we assume that the coupling is symmetric and we demonstrate that we can recover the strength of their coupling, K. In a second scenario, we do not assume symmetric coupling, and we recover the values of two coupling parameters, K 12 and K 21 (which should be identical). Then, we show that for 28 randomly coupled Rössler-like oscillators, the UKF technique can reconstruct the network architecture, inferring with high precision the presence or absence of links between pairs of oscillators. We investigate the quality of the network reconstruction across a wide range of values of coupling strength and show that we can accurately infer the network topology and the strength of the couplings, even when the oscillators are close to synchronization. This paper is organized as follows. The UKF algorithm and the model of the chaotic Rössler-like oscillators used in this study are described in section "Methods"; in section "Results" we perform different demonstrations using the algorithm to prove that we are able to characterize the topology of networks using noisy data. In section "Discussion" we present our conclusions and discuss the advantages and limitations of this technique.

Results
Uncoupled dynamics: estimation of the unobserved oscillators' variables. We begin by demonstrating a simple application of the UKF technique. We consider a synthetic data set corresponding to an isolated, uncoupled oscillator. We set N = 1 in Eqs. (7a-c) and solve the equations as described in "Methods" using the parameters listed in Table 1. Both, the numerically generated noisy values of the variables x, y, and z and the filtered time series obtained as described below are shown in Fig. 1. We remark that the UKF algorithm assumes that the input data are noisy measurements of an stochastic dynamical system, i.e., the input data is generated by a dynamically noisy system and it is contaminated by observational noise.
As expected, the UKF filters out the noise, recovering the smooth dynamics of the oscillator. This can be done even with high levels of observation noise (see Eq. 4). Therefore, the reconstruction of the evolution of the  www.nature.com/scientificreports/ three variables defining the state of the system was performed using only the noisy time series of the y variable. The algorithm may also be applied to observations of the other variables, x or z (i.e., from noisy observations of x we can reconstruct the evolution of the unobserved y and z). This confirms that the algorithm can infer the full state of the system from partial information of its dynamics. This feature, which is at the heart of the KF technique, introduces high versatility to the method presented here for systems for which the observations are only partial and noisy.
Minimal motif: two coupled oscillators. We now expand our analysis to a minimal motif formed by a system composed of two coupled oscillators. We set N = 2 and simulate Eqs. (7a-c) using a value of the coupling strength K that we will try to recover, from the analysis of the time series of the variables y 1 and y 2 of the two oscillators.
To estimate K, we extend the UKF state as shown in Eq. (5). We generated synthetic y 1 and y 2 measurement data time series for different coupling strengths K. In Fig. 2a, for example, we plot the time evolution of variables y 1 and y 2 when K = 0.1 . It can be seen that the oscillators are significantly coupled, but not fully synchronized. From these two time-traces, we can estimate both, the state of the two oscillators and K.
In order to apply the UKF, we prescribe 33 the model covariance 7 × 7 matrix, Q ω,e , (see Eq. 6) to be where Q ω is a covariance matrix estimated using a realization of a Gaussian process with σ ω = 0.02 , and the 2 × 2 measurement covariance Q ν is estimated using a realization of a Gaussian process with σ ν = 0.5 (see Table 1).
To run the UKF algorithm, a priori, we do know the exact value of the covariances. However, Q ω covariances created from σ ω ranging from 0.001 to 5 and Q ν covariances created from σ ν ranging from 0.001 to 1 have been tested. In all cases the algorithm converges to the correct value of K, showing its robustness to initial guesses.
(1) Q ω,e = Q ω 0 0 0  www.nature.com/scientificreports/ In Fig. 2b, the evolution of the estimate of K is shown. In this panel, we see that the real coupling strength (0.1) is obtained after a few iterations of the algorithm. Figure 2c shows the time evolution of variables y 1 and y 2 when K = 0.7 . In this case, the evolution of both elements is partially synchronized. The inference of K is shown in Fig. 2d. The parameter is accurately estimated (blue line) with less than 0.1% relative error, and the estimation is within one standard deviation (within the blue lighter lines) before 1/3 of the time series has been analysed. For the case of full synchronisation shown in Fig. 2e, we see that the algorithm, after a transient, also converges to the correct value (see Fig. 2f). The ability of the UKF algorithm to infer the correct value of the parameter despite the high level of synchronization of the two oscillators is a remarkable property worth noting.
This result confirms the ability of the UKF to recover correctly a single value of K. However, one may argue, as results shown in Fig. 2 seem to suggest, whether different levels of synchronization might make it more challenging to correctly estimate the coupling (specifically, from the noisy and partial measurement of the states of the oscillators). Therefore, to test the robustness of the UKF method, we created synthetic data sets for different values of K ranging from 0 to 1. The evolution of the estimations of the parameter K for these data sets is shown in Fig. 3. As it can be seen, the UKF is able to give a good estimate across the whole range of K values. We speculate that, if the length of the data-set used is increased, the result of the estimate will converge to the neighborhood a fixed value. One can also notice that the convergence takes more time as K increases. In all cases, the value of K = 0.5 was selected as the initial estimate because this is the value which is in the middle of the parameter range K = 0 − 1 . Since the initial guess for K is always the same, the reason for the increasingly long transient is not an increasing distance to the real value of K, thus we wonder that the longer transient may be related to the synchronization level of the signals. The parameter convergence depends on this initial estimate, however, the impact in our system of the initial estimate is minimal and does not significantly affect the parameter estimation performance or output. Figure 4 shows boxplot distributions of the UKF estimated K values versus the true K values used to generate the synthetic time series for y 1 and y 2 . To generate the figure, we considered the last 90,000 of the 100,000 data points obtained by the UK filter process. It is worth noting that the estimation for K = 0.6 is significantly away from the optimal result, while the estimates for the other K values are reasonably close to the target values (see also Fig. 3). In general, we can observe a worse agreement in the case in which the coupling is relatively high, where the oscillators start synchronizing.
Coupling directionality. So far, we assumed that the oscillators are symmetrically coupled (our prior knowledge is that the adjacency matrix is symmetric). We now relax this assumption by defining two coupling constants: K 12 = K · A 12 and K 21 = K · A 21 . By doing so, we check that the method can detect the directionality of the coupling.
We see in Fig. 5 that the UKF algorithm provides reliable estimates of K 12 and K 21 , but the estimates are more scattered, in comparison of the estimation shown in Fig. 4. Even though the UKF algorithm tends to slightly underestimate (overestimate) the true K 12 ( K 21 ) values, it provides good estimates for a wide range of Ks. One may note that the whiskers cover the diagonal dashed line in most cases, meaning that the estimates of the parameters are statistically valid. While we don't know why the estimates of the parameters K 12 and K 21 are more www.nature.com/scientificreports/ dispersed, we speculate that it may be due to the fact that we are estimating two parameters ( K 12 and K 21 ) instead of one (K), which possibly requires longer datasets for the parameters to converge, and/or the algorithm may be more sensitive to small fluctuations. Further tests are needed to determine how the sensitivity of the algorithm depends on the number of estimated parameters, the amount of data, the computational power required, and the quality of the results obtained.  www.nature.com/scientificreports/ Network of 28 randomly coupled Rössler-like chaotic oscillators. Lastly, we tackled the analysis of a network of N = 28 oscillators. We simulate Eqs. (7a-c) using a known value of the coupling strength K and a known adjacency matrix A ij (Network 1 in ref. 41 ). In this way, we obtain a synthetic data set with 28-time series ( y i with i = 1 . . . 28 ) which we use to first try to recover K (assuming that the adjacency matrix A ij is known) and then, we try to recover all the 378 ( N(N − 1)/2 ) values of K · A ij , assuming that A ii = 0 and A ij = A ji . Figure 6 presents the results obtained when inferring the value of K when we know the network topology, A ij . As it can be seen, the UKF quickly recovers the correct value of the coupling parameter. The broadest discrepancies are, again, for high values of K. Notice, also, that in these cases, the duration of the time traces are shorter than in Fig. 3. The reason for this is that, here, the level of noise is smaller (see Table 1).
The results obtained when estimating the 378 coupling parameters with K = 1 are displayed in Fig. 7a. It is clear that the algorithm can distinguish two different populations of K ij : one placed around 1 (the links) and another one around 0. Notice that, in this case, longer time traces are needed to obtain full convergence of all the estimated parameters. The reason for this is two-fold: the large number of parameters estimated, and the fact that  Table 1. In this case, the observation noise is smaller (which makes a faster convergence of the estimations) and the duration of the time traces is shorter than in Fig. 3. The solid black lines indicate the values of the K parameters to be estimated during the filtering. (b) Comparison of the estimated and synthetic K using boxplots. For each boxplot, the red line is the median, the upper and lower whiskers of the bars represent the maximum and minimum limits, and the red points refer to outliers. In these simulations, we have used fully stabilized parameters so they show almost perfect estimated values.  Table 1). The coupling strengths converge to either 0 or 1 representing not connection and connection, respectively, between nodes i and j. (b) The distance between the inferred adjacency matrix and the real one, D(K, K UKF ) , converges to zero. In the inset of panel (b) we show the logarithmic evolution of the distance, D, showing that all the estimated coupling terms converge to the correct value. www.nature.com/scientificreports/ we use K = 1 , which we have seen is more computationally demanding. In fact, the computational time needed here is about 5000 s contrasting with the few tens of seconds required for the previous results. The main reason for this change is the higher number of parameters to be estimated. It should be noted that no estimated parameter changes from one group of values to the other (from close to zero and close to one). However, they essentially approach the asymptotic value more or less monotonously. In other words, if we fix a threshold, say, equal to 0.5, promptly we can classify the parameter values into two groups. This allows us to obtain the adjacency matrix, A ij , without the burden of long simulations.
To quantify the quality of the network reconstruction we computed the Euclidean distance between the real and reconstructed matrices: The results are presented in Fig. 7b, where we see that D quickly tends to 0. By forcing the estimated parameters which are, say, above 0.5 to one and forcing the rest of the estimated parameters to zero, we perfectly reconstruct the real adjacency matrix.

Discussion and conclusions
To summarize, we have shown that the connectivity of a network of coupled oscillators can be inferred by using the Kalman filter technique. In order to do this, a good knowledge of the equations that describe the system is required. However, the level of confidence in this knowledge can be tuned by adjusting the matrix Q ω . From the experimental side, observed data is required (e.g., y i (t) for i = 1, . . . , N ) and a measurement function that relates the system variable(s) with the observed data is also needed. Again, the level of confidence in this knowledge can be tuned by adjusting the matrix Q ν .
Specifically, we have shown that the Unscented Kalman Filter (UKF) allows us to infer the connectivity by analyzing the dynamical evolution of only one of the three variables that define the phase space of each individual oscillator. After considering the state estimation of an isolated oscillator, we expanded the study to the minimal motif of two coupled oscillators. Using this approach we show that both the coupling strength and the directionality of connections can be recovered, even if the oscillators' signals are noisy. For this noisy signals, in order to obtain a better estimate of K, it is necessary to run the algorithm for much longer times, which is computationally demanding (especially for large values of N). We remark that the network topology can be inferred even when the oscillators are nearly synchronized.
We further tested the method by considering a random network of 28 oscillators. Again, we were able to obtain the coupling strength for the full range of values going from completely uncoupled, K = 0 , to fully coupled, K = 1 , with highly synchronized dynamics. In both situations, N = 2 and N = 28 , it is observed that the inference of K requires more data values for higher synchronization levels. However, if enough data is available, the correct estimation is obtained. As a final proof of the good performance of the method, the correct adjacency matrix was obtained for the network of N = 28 oscillators for K = 1. www.nature.com/scientificreports/ While the KF technique has important advantages, it also has some limitations. A main advantage is that Kalman filtering requires only partial observation of the state of the system to make predictions. This is an important feature of the method because the capability of inferring variables which are not accessible experimentally is remarkable in cases where experimental observations are limited. A drawback of the methodology is that it requires sufficiently long time traces for the convergence of the estimated values of the parameters. Another important drawback is that it requires a good knowledge of the system, i.e., a good model (with unknown parameters) needs to be used in the prediction step and a good measurement function linking the model to experimental observations is also needed.
As explained in the Introduction, we have used the implementation of Rössler oscillators given in Eqs. (7a-c) because it corresponds to an electronic implementation from which large experimental data sets are available 39,41 . Our objective is to be able to recover the experimental couplings. However, to do so, more information of the system is required. We have tried to infer connectivities from the experimental data sets presented in 39,41 ; however, the application of the UKF algorithm to simulations using Eqs. (7a-c) with the parameters indicated in 39,41 failed to provide a good agreement between the simulated time series and the experimentally measured ones. This may be due to typos in the parameter values and/or it might be due to the filtering of the signals performed by the data acquisition system. Improved knowledge of the experimental setup and data-acquisition methods is required for obtaining, using our method, a good reconstruction of the network.
It will be interesting, for future work, to apply the KF technique to uncover the connectivity of a network that evolves in time, where, for example, the strength of the connections varies. It will also be appealing to study the suitability of the KF technique when the observed data is heavily contaminated by noise. The results obtained here are promising and justify trying to use this approach in different scenarios.

Methods
Kalman filter. The KF, in its different implementations, is a standard tool used to estimate the unknown state of a system using a sequence of discrete-time observations and a model of the system evolution. Additionally, the filter can be used to infer a parameter value 33 by considering it as an additional system variable with the trivial dynamic evolution (being constant).
Originally, the filter was created for linear dynamical systems assuming additive Gaussian noise processes adding uncertainty to the state evolution and the observations. In this way, the algorithm makes recursive optimal estimates of the state of the system and its parameters 30 . The filter was soon extended to be able to deal with nonlinear systems 42 . Here, we use the Unscented Kalman Filter 43,44 implemented using the package unscented-KalmanFilter from Matlab.
The assumption underlying the UKF is that the system evolves as a noisy first-order Markov process where ū k corresponds to the state at the discrete-time k, ā is the vectorial function that maps the state at time k to the state at the next time, k + 1 , and ω k is the dynamic (or process) noise vector. The system is measured at each time k applying a measurement function b to the state ū k , perturbed by the measurement noise ν k : Both the dynamic noise ω k and the measurement noise ν k are additive zero-mean Gaussian processes with covariance Q ω and Q ν , respectively.
The algorithm works as follows. The estimate of the hidden state of the system u k at a given time is obtained in two steps. First, a new prediction is done, using ā and the previous estimate of the state. The information carried by the system observations, w k+1 , is then integrated with the prediction. The result is a prediction-correction algorithm that estimates the state ū k+1 . The algorithm also returns an estimate of the state covariance matrix, which provides the uncertainty on the estimation of ū k+1 .
The uncertainty of the model prediction and of the experimental information is represented by the covariances of the stochastic terms ω i and ν i . The filter is run using initial guesses of these two terms, so, by choosing them appropriately, the uncertainty of the model and of the experimental measurements is adjusted.
In the case in which q models parameters, (p 1 , . . . , p q ) have to be estimated, the hidden state is extended to incorporate them as additional variables: The trivial evolution of the parameters is governed by ṗ i = 0 , without process noise, resulting in an extended covariance matrix 33 : (5) u e = u 1 , . . . , u n , p 1 , . . . , p q www.nature.com/scientificreports/ This matrix is the extended form presented in Eq. (1), but now for q unknown parameters. The terms Q ω ij correspond to the elements of the covariance matrix obtained from the noise term ω k in Eq. (3) for the n dynamical variables of the system. We estimate these terms using a realization of the dynamical noise shown in Eq. (3). The rest of the terms of the matrix are zeros, corresponding to the covariance of the q constant parameters to be adjusted by the algorithm.
Models and data sets. In this work, we simulated networks composed of N Rössler-like oscillators described by Eqs. (1)-(4) in Sevilla-Escoboza et al 39,41 . The equations describing the system are the following: where We used the same values for all the parameters as in 39 , except for α 3 = 10,000 45 . Details about the value of the parameters can be found in Table 1. Experimental time series of the variable y i ( i = 1, . . . , N ) for each oscillator of the network are freely available for a wide range of coupling strengths and topologies 39,41 . Here, as a proof-ofconcept, we applied an UKF to simulated time series of the variables y i to infer the state of the oscillators and their connectivity under the assumption that we have a precise knowledge of the governing equations of the individual elements of the system. In the future, using these results as a starting point, we will try to reconstruct the networks from the available experimental data. However, this requires incorporating the manipulations performed during the data acquisition into the measurement function of the UKF (Eq. 4), since the available experimental data is not faithfully represented only by the raw solution of the Eqs. (1)-(4) in Ref. 39 .
To use the UKF, we transform this set of equations to the form defined in Eq. (3) which includes the presence of noise in the system ( ω i x (t k ), ω i y (t k ), ω i z (t k ) ). From this respect, the function ā is a Runge-Kutta step of fourth-order and a time-step dt defined in Table 1, while ω are Gaussian stochastic processes having zero mean and standard deviation σ ω (see Table 1). From realizations of these noise terms we estimate the covariance Q ω used in the filter algorithm.
Finally, to mimic a measurement process, we perturbed the y i variable with a Gaussian stochastic term as defined in Eq. (4) having zero mean and standard deviation σ ν (see Table 1). We decided to measure only variable y for two reasons: (i) in the experiments in 39 , y was the recorded variable for each oscillator, (ii) the coupling between the oscillators is reached by coupling their y i . This choice corresponds to b (x 1 , y 1 , z 1 , . . . , x N , y N , z N ) = (y 1 , y 2 , . . . , y N ) in Eq. (4).
In Eqs. (7a-c), the oscillators are mutually coupled with identical coupling strength, K. The coupling topology, for the different networks considered, is described by the adjacency matrix A defined as follows: A ij = A ji = 1 if oscillators i and j are connected, and A ij = A ji = 0 if they are not. There are no feedback loops, so we take A ii = 0 . In the case of N = 2 oscillators, to recover the coupling strength, we used the UKF to infer the value of one parameter only, K, and to reconstruct the directionality of the coupling, we used the UKF to infer the values of two parameters: K 12 = A 12 ⋅ K and K 21 = A 21 ⋅ K . For the analysis of N = 28 oscillators (in this case parameters are the same as in 41 ), we considered two different scenarios. First, we assumed that we know the adjacency matrix A, and we tried to infer the common coupling strength K, which is the same for all connections. In the second situation, we assumed bidirectionality and we tried to infer the N(N − 1)/2 = 378 unknown elements of K ⋅ A ij . (7a)