Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2010 October 7.
Published in final edited form as:
Phys Rev E Stat Nonlin Soft Matter Phys. 2009 May; 79(5 Pt 1): 051909.
Published online 2009 May 13. doi:  10.1103/PhysRevE.79.051909
PMCID: PMC2951269

Data assimilation for heterogeneous networks: The consensus set


Data assimilation in dynamical networks is intrinsically challenging. A method is introduced for the tracking of heterogeneous networks of oscillators or excitable cells in a nonstationary environment, using a homogeneous model network to expedite the accurate reconstruction of parameters and unobserved variables. An implementation using ensemble Kalman filtering to track the states of the heterogeneous network is demonstrated on simulated data and applied to a mammalian brain network experiment. The approach has broad applicability for the prediction and control of biological and physical networks.

Data assimilation has grown in importance for the analysis of experiments in physics, in part due to its great success in data-rich, high-dimensional problems. In modern numerical weather prediction, for example, it is common to use advanced filtering and prediction techniques for state estimation and tracking in cases where the experimental data are complicated and high dimensional [1].

Dynamical networks have become increasingly visible in physical and biological applications. Coupled grids of varying connectivity, often with simple dynamics at each grid node, can lead to interesting nonlinear phenomena, as seen in Josephson-junction arrays [2], neuronal networks [3], and a widening array of examples in systems biology [4]. The network structure has been found to be a decisive factor in the resulting dynamics, and is an area of intensive contemporary research [5].

The decentralized nature of networks presents increased difficulty for data assimilation. We propose a technique, called the consensus set method, that makes tracking and prediction of dynamics possible in a network where the equations of motion are unknown or only partly known, and when the system is only partially observed. After demonstrating the method computationally in heterogeneous model networks of excitable cells, we show the ability of the method to track experimental data from a mammalian brain network preparation for which no first-principles model exists.

In fact, the main targets of the technique are networks whose nodal dynamics is poorly described. In such cases, data assimilation may be attempted using models that only qualitatively approximate observed dynamics. Although model parameters and state information determined in this process may have little or no physical significance, the goal of the consensus set method is to capture and track global dynamics despite the lack of system knowledge. In particular, in the following we focus on tracking spiral wave and aperiodic behaviors in networks of excitable cells for prediction and control purposes.

The success of data assimilation is closely linked to the mathematical model for the underlying dynamics of the experiment. State-of-the-art methods such as ensemble Kalman filtering (EnKF) normally rely on good models to track the system and update the current system state (called the background in numerical weather prediction) [6]. The EnKF [711] is routinely used to complete partial state information. More recently, it has been extended to an algorithm for fitting a limited number of parameters as well as reconstructing the system state, in the presence of significant noise levels, sometimes called the dual-estimation problem [12,13]. The difficulty of fitting increases with the number of parameters needed, due to the rapid increase in data requirements.

Dynamical network models provide a demanding example of data requirements. It is typical for the number of model parameters to proliferate to an unmanageable number, since each element of the network may follow a slightly different local nonlinear dynamical model. Such a situation may overwhelm even the most robust parameter-fitting algorithms. We differentiate between two different types of network heterogeneities, calling a network weakly heterogeneous if the individual nodes of the network can be modeled by similar systems with differing parameter settings, and strongly heterogeneous if no such simplification exists. We will show the consensus set method, implemented below in the context of the unscented EnKF [13], to be useful for tracking and parameter identification in weakly heterogeneous networks undergoing complicated and nonstationary nonlinear dynamics, and in some cases for strongly heterogeneous networks.

First consider the weakly heterogeneous case. Begin with a general dynamical model with k parameters. Assume that several such versions of the model, with different parameter settings, are part of a network, and assume that the resulting dynamics is partially observable. In a moderately large network undergoing nonlinear dynamics, fitting all parameter settings in the network, along with reconstructing the full dynamical state from partial observations, will be infeasible. The idea of the consensus method is to track the weakly heterogeneous network with a homogeneous network of identical systems with a single set of k parameters that are fitted by the EnKF. The tracking will be inexact, because a heterogeneous network is being tracked by a homogeneous network. However, under the assumption that the heterogeneity is sufficiently weak, we show the ability of the EnKF to determine a consensus set of parameters that enables the homogeneous network to follow the heterogeneous network closely enough for many prediction and control purposes. The goal is not mathematical exactness, but an approximation that is close enough to be useful in an experimental scenario.

We begin by demonstrating the consensus set method in networks of excitable cells. The Wilson-Cowan equations [14] are designed to model the time evolution of cortical neurons with excitation u and recovery variable a. This system of two differential equations, with a single nonlinearity, supports sustained oscillations. A network of Wilson-Cowan oscillators on an N × N rectangular grid can generate a wide variety of interesting behaviors, including plane and ring waves, spiral waves, and chaotic aperiodic dynamics [15]. The equations at the (i, j) grid point are


where the sum is taken over all other nodes i′j′ in the network, d denotes the distance between grid points ij and i′j′, and H denotes the Heaviside function or a smoothed version if a continuous equation is desired [16]. The u variable represents instantaneous excitation, and the a variable represents recovery or adaptation. Biologically meaningful parameters ordinarily require cij1, cij2, cij6<0 and cij3, cij5>0.

Each model neuron in the grid follows Wilson-Cowan dynamics with input determined by superthreshold activity of its neighbors. The input is zero unless the variable ui′j′ is greater than the threshold cij4, in which case the neighbor contributes to the input according to its distance in the grid. Initial conditions are used that lead to a sustained nonequilibrium dynamical attractor on the grid. Observing only the u variable with applied noise at the grid points, we attempt to carry out approximate tracking of all u and a variables. To accomplish this, the EnKF will fit unknown parameters in the homogeneous Wilson-Cowan network model.

The Kalman filter uses known linear dynamical equations and observation functions along with observed data to continuously update a Gaussian approximation for the state and its uncertainty. The unscented EnKF [13] accomplishes these goals in the case of nonlinear dynamical equations. At each integration step, system states that are consistent with the current state uncertainty, called sigma points, are chosen. The EnKF consists of integrating the system from the sigma points, estimating mean state values, and then updating the covariance matrix that approximates the state uncertainty. The Kalman gain matrix updates the new most likely state of the system. As pointed out in recent literature [10,13] the ensemble Kalman filter can also be used to fit system parameters from data, introducing the unknown parameters as extra state variables with trivial dynamics. The EnKF with random initial conditions for the parameters will converge to the correct parameters, or in the case of slowly varying parameters, can track them along with the state variables. In recent work [17], we have shown the ability to track a parameter and spatiotemporal variables in homogeneous networks using EnKF in the presence of significant noise.

In the application of the unscented EnKF to an N × N network, the state is the dimension 2N2+m vector composed of the N2 states uij, the N2 states aij, and the m consensus parameters. The ensemble size is twice the dimension of the state vector, and we used an initial ensemble spread of 0.001 times the standard deviation of the variables aij (results were not sensitive to this choice). The one-step dynamics is the time-T map of the Wilson-Cowan system for some fixed time step T, together with the identity function in the m parameters.

As a first example, we consider a weakly heterogeneous 16 × 16 network of Wilson-Cowan neurons, where the parameters cijk are set randomly around the means ck, where c1=−0.12, c2=−0.64, c3=0.22, c4=0.4, c5=1, and c6=−0.056. A 10% heterogeneity is used, meaning that for each k=1, …, 6, the 256 coefficients cijk across the network have standard deviation equal to 10% of the mean ck. Most random choices of the parameters near these values support chaotic dynamics.

Figure 1 shows the result of solving for consensus parameters by the EnKF for this network, when the noisy versions of both uij and aij are observed. The observation noise added is 0.4 standard deviations of the clean signal, or 40% noise. Starting from random values, the parameters evolve over 15 s to stable values that roughly approximate the mean values of the cijk. From the experimental point of view, this is an uninteresting case, since we would like to observe only one voltage uij at each grid point ij, and consider aij as a hidden variable to be reconstructed, along with the parameters cijk. To carry out assimilation using only the observed voltages uij as input, we must address a hidden symmetry of model (1).

FIG. 1
The six parameters estimated by the EnKF converge over time to consensus values. The network's mean parameter values are denoted by closed circles. Here both the voltage and recovery variables u and a are observed.

Dual-estimation calculations often call for special treatment, because of the interaction of unknown variables and unknown parameters in the model. We next illustrate this problem in the case of the Wilson-Cowan system, and show a general approach that can be applied to directly handle experimental data. Consider the task of determining a(t) and c1, …, c6 in Eq. (1) from the knowledge of u(t) alone (we drop the ij subscripts for simplicity). Any experimental data u that are consistent with a(t), c1, c2, c3, c4, c5, c6 are consistent with a(t) / γ, c1, γc2, c3, c4, c5 / γ, c6 for arbitrary nonzero γ, meaning that the dual-estimation problem is underdetermined. Application of the EnKF directly to this problem results in parameter drift along the manifold of consistent solutions.

To achieve data assimilation in this context, we must force γ to assume a unique value. For example, if c5 is fixed in the EnKF to the same value as used in generating the input data, then γ=1 is forced, and the EnKF will succeed in fitting the remaining five parameters. Figure 2 shows the result of fixing c5=1 within the EnKF in the consensus set context. The input data are the same as in Fig. 1 except that only uij (with added noise) is observed by the EnKF. The five remaining consensus parameters are approximated well, similarly as in Fig. 1, in addition to the recovery variables aij. The variables at a typical ij grid point are shown in Fig. 2, verifying the ability of the method to successfully reconstruct aij.

FIG. 2
(Color online) The noisy uij input data and the output aij, shown as the black dashed line, from the consensus method at a typical grid point. The black solid line is the actual value of aij, for comparison. The red (light gray) dashed line is the filtered ...

In an experimental situation with poor prior knowledge of model dynamics, a correct value for c5 is too much to assume. We repeat the exercise with a 10% heterogeneous network after changing the mean of parameter cij5 to 0.5, leaving the other parameters as above. Since the EnKF has c5 fixed at 1, when presented with the u data it will reconstruct the symmetry with γ=0.5. As Fig. 3 shows, the result is the same output as in Fig. 2 except for aij replaced by aij / γ=2aij, and c2 replaced by γc2=−0.32, as expected from the above reasoning. This demonstrates the ability of the consensus set method to reconstruct viable hidden variables together with unknown parameters even in cases of significant model mismatch.

FIG. 3
(Color online) (a) The noisy uij data and the output aij from the consensus method (black dashed line) at one of the grid points. The black solid line is the actual value of aij. The reconstructed aij is magnified by 1 / γ=2. (b) Consensus parameters ...

Finally, we are prepared to apply the consensus set method to an experiment first described in [18]. Each element in the model will correspond to the volume of neural tissue imaged with a photodetector after the tissue was saturated with voltage-sensitive dye and illuminated. The signal from each photodetector thus represents the local mean field of transmembrane voltage from the ensemble of cells imaged, and is analogous to u in Eq. (1). Inhibition was blocked (with bicuculline), and our modeling reflects the lack of synaptic inhibition (a is a recovery variable, reflecting phenomena such as inactivation of sodium channels and potassium membrane repolarization) [16].

A 10 × 10 grid of voltage-sensitive dye readings were observed as a multidimensional time series. With c5=3 held fixed to avoid the symmetry, the consensus parameters converge (in a manner similar to that in Figs. 1 and and3)3) to c1=−1.08, c2=−2.82, c3=2.01, c4=−0.51, and c6=−0.20. Note that parameters c3 and c5 are positive, as anticipated in biologically reasonable models. Figure 4 shows the recovery variable aij created by the EnKF that is sufficiently compatible with the observed signal uij to allow tracking.

FIG. 4
(Color online) (a) Snapshot of the transmembrane voltage optical intensities uij from the neural tissue. (b) The experimental uij plotted as blue (dark gray) points, the filtered uij plotted as red (light gray) points, and the output aij from the consensus ...

In a larger sense, reconstructing equivalent dynamics may be more relevant and feasible goal than validating explicit model variables in situations where model inadequacy is severe [19]. Along similar lines, there has been recent work [20] emphasizing the role of dynamical synchronization over model fitting in data assimilation. The method outlined here can be used in conjunction with localized versions of EnKF [21] designed to address system spatial complexity as well as exploit parallel computation.

We view the consensus set method as a crucial step on the path to the eventual goal of prediction and control of complex dynamical networks, in particular when the individual network nodes are heterogeneous, and in the presence of model uncertainty. Some replacement for parameter fitting is essential for robust data assimilation in the face of model error, and this method, while not completely eliminating the problem of model error, may be a near-optimal choice in such a scenario.


We thank Richard Murray for helpful discussions. This research was supported by NSF Grant No. DMS-0811096, and NIH Grants No. R01MH50006 and No. K02MH01493.


1. Kalnay E. Atmospheric Modeling, Data Assimilation, and Predictability. Cambridge University Press; Cambridge, England: 2003.
2. Hadley P, Beasley MR, Wiesenfeld K. Phys Rev B. 1988;38:8712. [PubMed]Wiesenfeld K. Physica B. 1996;222:315.
3. Abbott LF, van Vreeswijk C. Phys Rev E. 1993;48:1483. [PubMed]Hopfield JJ, Herz AVM. Proc Natl Acad Sci USA. 1995;92:6655. [PubMed]
4. Jeong H, Tombor B, Albert R, Oltvai ZN, Barabasi AL. Nature (London) 2000;407:651. [PubMed]Jeong H, Mason SP, Barabasi AL, Oltvai ZN. ibid. 2001;411:41. [PubMed]
5. Newman MEJ. Phys Rev Lett. 2002;89:208701. [PubMed]Pecora LM, Carroll TL. ibid. 1998;80:2109.Barahona M, Pecora LM. ibid. 2002;89:054101. [PubMed]Nishikawa T, Motter AE, Lai YC, Hoppensteadt FC. ibid. 2003;91:014101. [PubMed]Donetti L, Hurtado PI, Munoz MA. ibid. 2005;95:188701. [PubMed]
6. Kalman R. ASME J Basic Eng. 1960;82:35.
7. Evensen G. J Geophys Res. 1994;99:10143.
8. Evensen G, Van Leeuwen P. Mon Weather Rev. 2000;128:1852.
9. Hunt B, Kalnay E, Kostelich E, Ott E, Patil DJ, Sauer T, Szunyogh I, Yorke JA, Zimin A. Tellus, Ser A. 2004;56:273.
10. Hansen JA, Penland C. Physica D. 2007;230:88.
11. Hunt BR, Kostelich EJ, Szunyogh I. Physica D. 2007;230:112.
12. Wan EA, Van der Merwe R, Nelson AT. Adv Neural Inf Process Syst. 2000;12:666.
13. Voss HU, Timmer J, Kurths J. Int J Bifurcation Chaos Appl Sci Eng. 2004;14:1905.
14. Wilson HR, Cowan JD. Biophys J. 1972;12:1. [PubMed]
15. Ermentrout GB. SIAM J Appl Math. 1984;44:80.Minai AA, Anand T. Phys Rev E. 1998;57:1559.Rabinovich MI, Huerta R, Varona P, Afraimovich VS. PLOS Comput Biol. 2008;4:e1000072. [PubMed]
16. Pinto D, Ermentrout GB. SIAM J Appl Math. 2001;62:206.
17. Schiff SJ, Sauer T. J Neural Eng. 2008;5:1. [PMC free article] [PubMed]
18. Huang X, Troy W, Yand Q, Ma H, Laing CR, Schiff SJ, Wu JY. J Neurosci. 2004;24:9897. [PubMed]Schiff SJ, Huang X, Wu JY. Phys Rev Lett. 2007;98:178102. [PubMed]
19. Baek SJ, Hunt B, Kalnay E, Ott E, Szunyogh I. Tellus, Ser A. 2006;58:293.
20. Duane GS, Tribbia JJ, Weiss JB. Nonlinear Processes Geophys. 2006;13:601.Yang SC, Baker D, Li H, Cordes K, Huff M, Nagpal G, Okereke E, Villafane J, Kalnay E, Duane G. J Atmos Sci. 2006;63:2340.
21. Ott E, Hunt BR, Szunyogh I, Zimin AV, Kostelich EJ, Corazza M, Kalnay E, Patil DJ, Yorke JA. Tellus, Ser A. 2004;56:415.