Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2951269

Formats

Article sections

Authors

Related links

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.051909PMCID: PMC2951269

NIHMSID: NIHMS133872

See other articles in PMC that cite the published article.

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 [7–11] 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

$$\begin{array}{l}{\dot{u}}_{ij}={c}_{ij1}{u}_{ij}+{c}_{ij2}{a}_{ij}+{c}_{ij3}\sum _{i\prime j\prime}{e}^{-d}H({u}_{i\prime j\prime}-{c}_{ij4}),\\ \phantom{\rule{5em}{0ex}}{\dot{a}}_{ij}={c}_{ij5}{u}_{ij}+{c}_{ij6}{a}_{ij},\end{array}$$

(1)

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 *c _{ij}*

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 *u _{i′j′}* is greater than the threshold

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 2*N*^{2}+*m* vector composed of the *N*^{2} states *u _{ij}*, the

As a first example, we consider a weakly heterogeneous 16 × 16 network of Wilson-Cowan neurons, where the parameters *c _{ijk}* are set randomly around the means

Figure 1 shows the result of solving for consensus parameters by the EnKF for this network, when the noisy versions of both *u _{ij}* and

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 *c*_{1}, …, *c*_{6} 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*)*, c*_{1}, *c*_{2}, *c*_{3}, *c*_{4}, *c*_{5}, *c*_{6} are consistent with *a*(*t*) / *γ*, *c*_{1}, *γc*_{2}, *c*_{3}, *c*_{4}, *c*_{5} / *γ*, *c*_{6} 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 *c*_{5} 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 *c*_{5}=1 within the EnKF in the consensus set context. The input data are the same as in Fig. 1 except that only *u _{ij}* (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

(Color online) The noisy *u*_{ij} input data and the output *a*_{ij}, shown as the black dashed line, from the consensus method at a typical grid point. The black solid line is the actual value of *a*_{ij}, 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 *c*_{5} is too much to assume. We repeat the exercise with a 10% heterogeneous network after changing the mean of parameter *c _{ij}*

(Color online) (a) The noisy *u*_{ij} data and the output *a*_{ij} from the consensus method (black dashed line) at one of the grid points. The black solid line is the actual value of *a*_{ij}. The reconstructed *a*_{ij} 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 *c*_{5}=3 held fixed to avoid the symmetry, the consensus parameters converge (in a manner similar to that in Figs. 1 and and3)3) to *c*_{1}=−1.08, *c*_{2}=−2.82, *c*_{3}=2.01, *c*_{4}=−0.51, and *c*_{6}=−0.20. Note that parameters *c*_{3} and *c*_{5} are positive, as anticipated in biologically reasonable models. Figure 4 shows the recovery variable *a _{ij}* created by the EnKF that is sufficiently compatible with the observed signal

(Color online) (a) Snapshot of the transmembrane voltage optical intensities *u*_{ij} from the neural tissue. (b) The experimental *u*_{ij} plotted as blue (dark gray) points, the filtered *u*_{ij} plotted as red (light gray) points, and the output *a*_{ij} 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |