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

**|**PLoS Comput Biol**|**v.6(5); 2010 May**|**PMC2865517

Formats

Article sections

- Abstract
- Author Summary
- Introduction
- Results
- Discussion
- Materials and Methods
- Supporting Information
- References

Authors

Related links

PLoS Comput Biol. 2010 May; 6(5): e1000776.

Published online 2010 May 6. doi: 10.1371/journal.pcbi.1000776

PMCID: PMC2865517

Karl J. Friston, Editor^{}

University College London, United Kingdom

* E-mail: ude.usp@minahg

Conceived and designed the experiments: GU SJS. Performed the experiments: GU. Analyzed the data: GU. Wrote the paper: GU SJS.

Received 2010 January 19; Accepted 2010 April 1.

Copyright Ullah, Schiff.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are properly credited.

This article has been cited by other articles in PMC.

Observability of a dynamical system requires an understanding of its state—the collective values of its variables. However, existing techniques are too limited to measure all but a small fraction of the physical variables and parameters of neuronal networks. We constructed models of the biophysical properties of neuronal membrane, synaptic, and microenvironment dynamics, and incorporated them into a model-based predictor-controller framework from modern control theory. We demonstrate that it is now possible to meaningfully estimate the dynamics of small neuronal networks using as few as a single measured variable. Specifically, we assimilate noisy membrane potential measurements from individual hippocampal neurons to reconstruct the dynamics of networks of these cells, their extracellular microenvironment, and the activities of different neuronal types during seizures. We use reconstruction to account for unmeasured parts of the neuronal system, relating micro-domain metabolic processes to cellular excitability, and validate the reconstruction of cellular dynamical interactions against actual measurements. Data assimilation, the fusing of measurement with computational models, has significant potential to improve the way we observe and understand brain dynamics.

To understand a complex system such as the weather or the brain, one needs an exhaustive detailing of the system variables and parameters. But such systems are vastly undersampled from existing technology. The alternative is to employ realistic computational models of the system dynamics to reconstruct the unobserved features. This model based state estimation is referred to as data assimilation. Modern robotics use data assimilation as the recursive predictive strategy that underlies the autonomous control performance of aerospace and terrestrial applications. We here adapt such data assimilation techniques to a computational model of the interplay of excitatory and inhibitory neurons during epileptic seizures. We show that incorporating lower scale metabolic models of potassium dynamics is essential for accuracy. We apply our strategy using data from simultaneous dual intracellular impalements of inhibitory and excitatory neurons. Our findings are, to our knowledge, the first validation of such data assimilation in neuronal dynamics.

A universal dilemma in understanding the brain is that it is complex, multiscale, nonlinear in space and time, and we never have more than partial experimental access to its dynamics. To better understand its function one not only needs to encompass the complexity and nonlinearity, but also estimate the unmeasured variables and parameters of brain dynamics. A parallel comparison can be drawn in weather forecasting [1], although atmospheric dynamics are arguably less complex and less nonlinear. Fortunately, the meteorological community has overcome some of these issues by using model based predictor-controller frameworks whose development derived from computational robotics requirements of aerospace programs in 1960s [2], [3]. A predictor-controller system employs a computational model to observe a dynamical system (e.g. weather), assimilate data through what may be relatively sparse sensors, and reconstruct and estimate the remainder of the unmeasured variables and parameters in light of available data. The result of future measured system dynamics is compared with the model predicted outcome, the expected errors within the model are updated and corrected, and the process repeats iteratively. For this recursive initial value problem to be meaningful one needs computational models of high fidelity to the dynamics of the natural systems, and explicit modeling of the uncertainties within the model and measurements [3]–[5].

The most prominent of the model based predictor-controller strategies is the Kalman filter (KF) [2]. In its original form, the KF solves a linear system. In situations of mild nonlinearity, the extended forms of the KF were used where the system equations could be linearized without losing too much of the qualitative nature of the system. Such linearization schemes are not suitable for neuronal systems with nonlinearities of the scale of action potential spike generation. With the advent of efficient nonlinear techniques in the 1990s such as the ensemble Kalman filter [6], [7] and the unscented Kalman filter (UKF) [8], [9], along with improved computational models for the dynamics of neuronal systems (incorporating synaptic inputs, cell types, and dynamic microenvironment) [10], the prospects for biophysically based ensemble filtering from neuronal systems are now strong. The general framework of the UKF differs from the extended KF in that it integrates the fundamental nonlinear models directly, along with iterating the error and noise expectations through these nonlinear equations. Instead of linearizing the system equations, UKF performs the prediction and update steps on an ensemble of potential system states. This ensemble gives a finite sampling representation of the probability distribution function of the system state [3], [11]–[15].

Our hypothesis is that seizures arise from a complex nonlinear interaction between specific excitatory and inhibitory neuronal sub-types [16]. The dynamics and excitability of such networks are further complicated by the fact that a variety of metabolic processes govern the excitability of those neuronal networks (such as potassium concentration () gradients and local oxygen availability), and these metabolic variables are not directly measurable using electrical potential measurements. Indeed, it is becoming increasingly apparent that electricity is not enough to describe a wide variety of neuronal phenomena. Several seizure prediction algorithms, based only on EEG signals, have achieved reasonable accuracy when applied to static time-series [17]–[19]. However, many techniques are hindered by high false positive rates, which render them unsuitable for clinical use. We presume that there are aspects of the dynamics of seizure onset and pre-seizure states that are not captured in current models when applied in real-time. In light of the dynamic nature of epilepsy, an approach that incorporates the time evolution of the underlying system for seizure prediction is required. As one cannot see much of an anticipatory signature in EEG dynamics prior to seizures, the same can be said of a variety of oscillatory transient phenomena in the nervous system ranging from up states [20], spinal cord burst firing [21], cortical oscillatory waves [22], in addition to animal [23] and human [24] epileptic seizures. All of these phenomena share the properties that they are episodic, oscillatory, and have apparent refractory periods following which small stimuli can both start and stop such events.

It has recently been shown that the interrelated dynamics of and sodium concentration () affect the excitability of neurons, help determine the occurrence of seizures, and affect the stability of persistent states of neuronal activity [10], [25]. Competition between intrinsic neuronal ion currents, sodium-potassium pumps, glia, and diffusion can produce slow and large-amplitude oscillations in ion concentrations similar to what is observed physiologically in seizures [26], [27].

Brain dynamics emerge from within a system of apparently unique complexity among the natural systems we observe. Even as multivariable sensing technology steadily improves, the near infinite dimensionality of the complex spatial extent of brain networks will require reconstruction through modeling. Since at present, our technical capabilities restrict us to only one or two variables at a restricted number of sites (such as voltage or calcium), computational models become the “lens” through which we must consider viewing all brain measurements [28]. In what follows, we will show the potential power of fusing physiological measurements with computational models. We will use reconstruction to account for unmeasured parts of the neuronal system, relating micro-domain metabolic processes to cellular excitability, and validating cellular dynamical reconstruction against actual measurements.

As a first example of assimilating neural data we used intracellular voltage data from a spiking pyramidal cell (PC) from the Cornu Ammonis region 1 (CA1) of rat hippocampus. Using only the noisy membrane potential measurement, , we employed modified Hodgkin-Huxley equations to reconstruct and track all of the gating variables of the ion channels: sodium channel activation and inactivation variables and , and potassium channel activation variable (Figure 1). Beginning with arbitrary initial conditions the root mean square (RMS) error between measured and estimated membrane potential changes with time (Figure 2). As is clear from the figure the RMS error converges to near zero within a few hundred milliseconds for the simulations shown in Figure 1. We also tracked the maximum conductance parameters of the ion channels (not shown).

Model inadequacy is an issue of intense research in the data assimilation community – no model does exactly what nature does. To deal with inadequate models, researchers in areas such as meteorology have developed various strategies to account for the inaccuracies in the models for weather forecasting [4], [5], [29]. In complex systems such as neuronal networks, the need to account for model inadequacy is critical. To demonstrate that UKF can track neuronal dynamics in the face of moderate inadequacy, we impaired our model by setting the sodium current rate constant instead of using the actual complex function of , (see equation (2) for the functional form of ), and tracked it as a parameter (Figure 3). That is, we deleted the relevant function for from the model and allowed UKF to update it as a parameter. The model with fixed is by itself unable to spike, but when it is allowed to float when voltage is assimilated through UKF using the data from hippocampal pyramidal cells (PCs), it is capable of tracking the dynamics of the cell reasonably well. The tracked by the filter is sufficiently close to its functional form values (within 25%) so that spiking dynamics can be reconstructed (Figure 3C and 3D). This occurs because Kalman filtering constantly estimates the trade off between model accuracy and measurements, expressed in the filter gain function [2], [3]. This is an excellent demonstration of the robustness of this framework. Looking at the estimated values of it also becomes clear that in fact should be assigned the functional form rather than a constant value.

Despite decades of effort neuroscientists lack a unifying dynamical principle for epilepsy. An incomplete knowledge of the neural interactions during seizures makes the quest for unifying principles especially difficult [30]. Here we show that UKF can be employed to track experimentally inaccessible neuronal dynamics during seizures. Specifically, we used UKF to assimilate data from pairs of simultaneously impaled pyramidal cells and oriens-lacunosum moleculare (OLM) interneurons (INs) in the CA1 area of the hippocampus [23]. We then used biophysical ionic models to estimate extra- and intracellular potassium, sodium, and calcium ion concentrations and various parameters controlling their dynamics during seizures (Figure 4). In Figure 4A we show an intracellular recording from a pyramidal cell during seizures, and plot the estimated extracellular potassium concentration () in Figure 4B. As is clear from the figure the extracellular potassium concentration oscillates as the cell goes into and out of seizures. The potassium concentration begins to rise as the cell enters seizures and peaks with the maximal firing frequency, followed by decreasing potassium concentration as the firing rate decreases and the seizure terminates. Higher makes the PC more excitable by raising the reversal potential for currents (equation 7). The increased reversal potential causes the cell to burst-fire spontaneously. Whether the increased causes the cells to seize or is the result of seizures has been an old question [31] whose resolution will likely take place from better understanding of the coupled dynamics. For present purposes, it is known that increased in experiments can support the generation, and increase the frequency and propagation velocity of seizures [32], [33]. Changes in the concentration of intracellular sodium ions, , are closely coupled with the changes of (Figure 4C). As shown in panels (4D–F) we reconstructed the parameters controlling the microenvironment of the cell. These parameters included the diffusion constant of in the extracellular space, buffering strength of glia, and concentration in the reservoir of the perfusing solution *in vitro* (or in the vasculature *in vivo*) during seizures. Note that the ionic concentration in the distant reservoir is different from the more rapid dynamics within the smaller connecting extracellular space near single cell where excitability is determined. We were also able to track other variables and parameters such as extracellular calcium concentration and ion channel conductances.

In Figure 5, we show an expanded view of a single cell response during a single seizure from Figure 4. Extracellular potassium concentration increases several fold above baseline values during seizures [31]. During a single seizure, starts rising from a baseline value of 3.0mM as the seizure begins and peaks at 7mM at the middle of the seizure (Figure 5). Interestingly the estimated by UKF matches very closely the measured seen *in vitro* studies [34].

Considering the slow time scale of seizure evolution (period of more than 100 seconds in our experiments), we test the importance of slow variables such as ion concentrations for seizure tracking. As shown in Figure 6, we found that including the dynamic intra- and extracellular ion concentrations in the model is necessary for accurate tracking of seizures. Using Hodgkin-Huxley type ionic currents with fixed intra- and extracellular ion concentration of and ions fails to track seizure dynamics in pyramidal cells (Figure 6C). We used physiologically normal concentrations of 4mM and 18mM for extracellular and intracellular respectively for these simulations. The conclusion remains the same when higher and are used. A similar tracking failure is found while tracking the dynamics of OLM interneurons during seizures (not shown). To further emphasize the importance of ion concentrations dynamics for tracking seizures we calculate the Akaike's information criterion (AIC) for the two models used in Figure 6, i.e. the model with and without ion concentration dynamics. AIC is a measure of the goodness of fit of a model and offers a measure of the information lost when a given model is used to describe experimental observations. Loosely speaking, it describes the tradeoff between precision and complexity of the model [35]. We used equation (29) for the AIC measure. The AIC measure for the model without ion concentration dynamics is . The model with ion concentration dynamics on the other hand has AIC value equal to , indicating the importance of ion concentration dynamics for tracking seizures.

Pyramidal cells and interneurons in the hippocampus reside in different layers with different cell densities. To investigate whether there exist significant differences in the microenvironment surrounding these two cell types we assimilated membrane potential data from OLM interneurons in the hippocampus and reconstructed and ion concentrations inside and outside the cells. As shown in Figure 7, both the baseline level and peak near the interneurons must be very high as compared to that seen for the pyramidal cells (cf. Figure 4B). This is an important prediction in light of the recently observed interplay between pyramidal cells and interneurons during *in vitro* seizures [23]; in these experiments pyramidal cells were silent when the interneurons were intensively firing. Following intense firing the interneurons entered a state of depolarization block simultaneously with the emergence of intense epileptiform firing in pyramidal cells. Such a novel pattern of interleaving neuronal activity is proposed to be a possible mechanism for the sudden drop in inhibition during seizures – it may be permissive of runaway excitatory activity. The mechanism leading to such interplay, specifically the reasons for differential firing patterns in pyramidal cells and interneurons are unknown. Our results here indicate the potential role of the neuronal microenvironment in producing such interplay. Our findings suggest that the buffering mechanism in the OLM layer is weaker as compared with the pyramidal layer, thus causing higher in the OLM layer. The higher surrounding the interneurons causes increased excitability of the cell by raising the reversal potential for currents (higher than the pyramidal cells, see equation 7). The higher reversal potential for currents causes the interneuron to spontaneously burst fire at higher frequency and eventually drives the interneuron to transition into depolarization block when firing is peaked. As the INs enter the depolarized state, the inhibitory synaptic input from the INs to the PCs drops substantially, releasing PCs to generate the intense excitatory activity of seizures (equation 8, Figure S3). The collapse of inhibition due to the entrance of INs into a depolarized state also helps explain the sudden decrease in inhibition at seizure onset in neocortex described by Trevelyan, et al. [36] as the loss of *inhibitory veto*. As shown in Figure S1, we also tracked the remaining variables for the INs.

Since the interaction of neurons determines network patterns of activity, it is within such interactions that we seek unifying principles for epilepsy. To demonstrate that the UKF framework can be utilized to study cellular interactions, we reconstructed the dynamics of one cell type by assimilating the measured data from another cell type in the network. In Figure 8 we only show the estimated membrane potentials, but we also reconstructed the remaining variables and parameters of both cells (Figures S2 and S3). We first assimilated the membrane potential of the PC to estimate the dynamics of the same cell and also the dynamics of a coupled IN (Figure 8A–D). Conversely, we estimate the dynamics of PC from the simultaneously measured membrane potential measurements of the IN (Figure 8D–F). As is evident from Figure 8 the filter framework is successful at reciprocally reconstructing and tracking the dynamics of these different cells within this network. In Figure S2, we show intracellular concentration and gating variables of and channels in PCs for simulation in Figure 8A–D. The variables modeling the synaptic inputs for both INs and PCs in Figure 8A–D are shown in Figure S3. As clear from Figure S3 (D), the variable (equation 8) reaches very high values when the INs lock into depolarization block, shutting off the inhibitory inputs from INs to PCs.

There has been intense interest in the neuroscience communities in bringing control-theoretical tools to bear on neuronal encoding and decoding problems [37]–[45]. In all of this work, statistical models (continuous or point process) were fit to data recorded from neurons, and these empirical models incorporated into applications. Our use of control theoretic tools is very different. We built computational models from the physiological properties of neurons and their networks, as well as the properties of ion metabolism, *without data fitting*. Using these fundamental models of the physics of neuronal systems, we fuse these models with data – *data assimilation* – in a manner commonly applied in meteorology [1], [46]–[50]. We are aware of a recent laboratory demonstration in fluid mechanics using a simplified model of fluid dynamics (Boussinesq equations) in a similar manner as we have performed here [51] (see also [14]). Other authors have also recently discussed the importance and power of going beyond statistical empirical models in neuronal systems, and simulations have begun to explore the feasibility of carrying this out [52]–[54]. To our knowledge, our findings are the first experimental validation that a fundamental biophysical model of part of the brain can be employed to assimilate incomplete data and accurately reconstruct its network dynamics.

Our conjecture is that the parallels with numerical meteorology are deep. By the turn of the 20th century, it was apparent that the lack of periodicities in weather limited forecasts based on previous state (autoregressive) statistical models, and that integrating the actual equations of motion of the atmosphere would be required. Infeasible initially, the turning point came when integrating such models gave ‘first approximations that bore a recognizable resemblance to the actual motions’ [55]. Furthermore, the use of simplified dynamical models that retained the most important of the physical dynamics was a critical development [1].

Our findings suggest that an analogous use of biophysical models of neuronal processes using the recursive predictive strategies employed in meteorological data assimilation is now feasible. We are presently exploring such application in frameworks for model-based data assimilation and control of Parkinson's disease [15]. Experiments are underway exploring the application for seizures in the intact brain, and the assimilation of cognitive rhythms. The potential for such techniques to improve our understanding of the dynamics of single cells and neuronal networks is substantial.

In conclusion, we have demonstrated the feasibility for data assimilation within neuronal networks using detailed biophysical models. In particular, we demonstrated that estimating the neuronal microenvironment and neuronal interactions can be performed by embedding our improving biophysical neuronal models within a model based state estimation framework. This approach can provide a more complete understanding of otherwise incompletely observed neuronal dynamics during normal and pathological brain function.

We used two-compartmental models for the pyramidal cells and interneurons: a cellular compartment and the surrounding extracellular microenvironment. The membrane potentials of both cells were modeled by Hodgkin-Huxley equations containing sodium, potassium, calcium-gated potassium (after-hyperpolarization), and leak currents. For the network model, the two cell types are coupled synaptically and through diffusion of potassium ions in the extracellular space. A schematic of the model is shown in Figure 9.

The membrane potential of the neurons is modeled with the following set of modified Hodgkin-Huxley equations [10],

(1)

where and represent gating variables for potassium, , and sodium, , currents. The leak current, , has three components: leak, , leak, , and chloride leak, . The after-hyperpolarization current is only included in the pyramidal cell to account for its frequency adaptation. The meaning and values of parameters used in the model are given in Table 1.

The rate equations for the gating variables are

(2)

The current equations were augmented with dynamic variables representing the intra- and extracellular ion concentrations (, , and ). These ion concentrations are affected by the neuron's intrinsic ionic currents as well as a sodium-potassium pump current. The glial buffering, diffusion between the nearest neighbor cells, and diffusion into the environment of the cell (bath solution in slice preparation and vasculature *in vivo*) also modulate the potassium concentration in the microenvironmental extracellular space. The ion concentrations inside and outside the cell are coupled to the membrane voltage equations via the Nernst equation [10], [13], [25]. Finally, PCs and INs are connected to each other through synaptic inputs and diffusion of extracellular potassium between the nearest neighbor neurons.

Given the potassium ion currents , activity of the pump exchanging and , , diffusion of potassium to the microenvironment, , and glial buffering, , the extracellular potassium dynamics, , can be represented in the model as (Figure 9).

(3)

where the − pump is modeled as a product of a sigmoidal functions, is the pump strength under normal conditions, and is the intracellular sodium concentration. Each sigmoidal term saturates for high values of internal sodium and external potassium respectively. More biophysically realistic models of pumps, such as those in [57] produce substantially similar results. in the diffusion equation is the potassium concentration in the nearby reservoir. Physiologically, this would correspond to either the bath solution in a slice preparation, or the vasculature in the intact brain (noting that is kept below the plasma level by trans-endothelial transport). Both active and passive uptake into glia is incorporated into a simplified single sigmoidal response function that depends on extracellular concentration with representing the maximum buffering strength. A similar but more physiological approach was used in [58]. Two factors allow the glia to provide a nearly insatiable buffer for the extracellular space. The first is the large size of the glial network. Second, the glial endfeet surround the pericapillary space, which, through interaction with arteriole walls, can effect blood flow; this in turn can increase the buffering capability of the glia [59]–[61].

We consider a spherical cell with a radius of . The diffusion coefficient of to the nearby reservoir , is obtained from Fick's law, /, where /sec is the diffusion constant in neocortex [62], and for brain reflects the average distance between capillaries [63]. The factor 0.3mM/coul in equation (3) converts ionic current to concentration rate of change and is calculated using / [10], where , and represent cell area, volume and Faraday constant respectively. is equal to 3mM in physiological conditions, and the intra- to extracellular volume ratio [64].

To complete the description of dynamics, we make the assumption that the flow of into the cell is compensated by flow of out of the cell to maintain bulk electroneutrality. Thus the internal potassium concentration () can be approximated by [10]

(4)

where 140mM and 18mM respectively are the normal resting concentrations of and inside the cell.

The intracellular and extracellular concentrations (, ) are also updated in the model as [10]

(5)

where 144mM is the normal resting extracellular concentration. in equations (3) and (5) are multiplied by factor 2 and 3 respectively due to the fact that the pump has an electrogenic 23 ratio.

The reversal potentials for , and are updated based on the instantaneous ion concentrations using the Nernst equations

(7)

Equation (7) binds the ion concentrations dynamics to the Hodgkin-Huxley equations (1, 2).

The pyramidal cells and OLM interneurons are coupled both synaptically and through extracellular diffusion as shown in Figure 9. The following synaptic currents are added to the membrane potential equations [25]

(8)

Where the superscripts and represent pyramidal cell and interneuron respectively. and is the membrane potential of the PCs and INs respectively. The variable takes into account the firing interplay between pyramidal cells and interneurons [25]. Ziburkus, et al. [23] observed during *in vitro* seizures that pyramidal cells were silent when the interneurons were burst firing, followed by high frequency firing in pyramidal cells when interneurons were locked into a depolarized state called depolarization block. The variable in equation (8) causes the synaptic input to drop to zero when the cells go to depolarization block. Various parameters used in equation (8) are: , , , , and . Synaptic strengths , are mimicked according to and AMPA inputs and values of 0.84 and 0.17, respectively, are used for the simulations. The variable gives the temporal evolution of the synaptic input from the pyramidal cell to the interneuron and is the synaptic input from the interneuron to the pyramidal cell. and evolve as

(9)

The parameters and are the time constants for the excitatory and inhibitory synapses respectively and .

In the case of coupled pyramidal cells and interneurons, the rate equation for is updated by adding the following lateral diffusion term ()

(10)

where is the separation between the interneurons and pyramidal cells.

To estimate and track the dynamics of the neuronal networks, we applied a nonlinear ensemble version of the Kalman filter, the unscented Kalman filter (UKF) [8], [9]. The UKF uses known nonlinear dynamical equations and observation functions along with noisy, partially observed data to continuously update a Gaussian approximation for the neuronal state and its uncertainty. At each integration step, perturbed system states that are consistent with the current state uncertainty, *sigma points*, are chosen. The UKF 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 based on the estimated measurements and the actual partially measured state. The estimated states (filtered states) are used to estimate the experimentally inaccessible parameters and variables by synchronizing the model equations to the estimated states. To estimate the system parameters from data, we introduced the unknown parameters as extra state variables with trivial dynamics. The UKF with random initial conditions for the parameters will converge to an optimal set of parameters, or in the case of varying parameters, will track them along with the state variables [11]–[13].

Given a function describing the dynamics of the system (equations 1–10 in our case), and an observation function contaminated by uncertainty characterized in the covariance matrix , for a -dimensional state vector with mean the UKF generates the
*sigma points*
, …, so that their sample mean and sample covariance are and . The *sigma points* are the rows of the matrix

(11)

The index on the left-hand side corresponds to the row taken from the matrix in the parenthesis on right-hand side. The square root sign denotes the matrix square root and indicates transpose of the matrix. Sigma points can be envisioned as sample points at the boundaries of a covariance ellipsoid. In what follows, superscript tilde ( ) represents the *a priori* values of variables and parameter, i.e. the values at a given time-step when observation up to time-step are available, while hat ( ) represents the *a posteriori* quantities, i.e. the values at time-step when observations up to time-step are available.

Applying one step of the dynamics to the *sigma points* and calling the results , and denoting the observations of the new states by , we define the means

(12)

where and are the *a priori* state and measurement estimates, respectively. Now define the *a priori* covariances

(13)

of the ensemble members. The Kalman filter estimates of the new state and uncertainty are given by the *a posteriori* quantities

(14)

and

(15)

where is the Kalman gain matrix and is the actual observation [3], [8], [9], [11]–[13]. Thus and are the updated estimated state and covariance for the next step. The *a posteriori* estimate of the observation is recovered by . Thus by augmenting the observed state variables with unobserved state variables and system parameters, UKF can estimate and track both unobserved variables and system parameters.

In our simulations, the state is the dimensional vector consisting of the variable values (equations 1–10) describing the dynamics of neurons and the parameter values to be tracked. The one-step dynamics function is the system of differential equations (equations 1–10). State vector for a single PC is given as

(16)

Where , , …. are the parameters that we want to track. For example, we tracked three parameters in Figure 4, replacing , , …. by , , and respectively in equation (16). For coupled PC and IN, the state vector included variables , , , , , , and for IN along with four synaptic variables, , , , and in order to represent the synaptic interactions between the two cells. The observation function returned the first component of the vector (membrane potential, ) at given time t. We observed the membrane potential of the cell and treated the rest of the variables as unobserved. For most of our simulations we used an integration time-step dt=0.01ms while the membrane potential of the neuron was measured each 0.1ms.

An iteration of the filter is performed in the following three steps (see [3], [8], [9] for more details).

*Initialization:* The filter is initialized as follows

(17)

where are the initial values of the state variables, and represent expectation.

*Prediction:* The following equations are used to propagate the state estimate and covariance from time-step (k−1) to k. First create a set of *sigma points* by applying equation (11) to system state equation (16)

(18)

The *sigma points* are transformed into vectors by using the nonlinear system of equations (1–10)

(19)

The average of vectors gives the *a priori* state estimate at time .

(20)

The *a priori* error covariance is given as

(21)

where represents the process noise.

*Measurement Update:* We implemented the measurement update as follows. Given the current guess for the mean, , and covariance, of , we choose new sigma points

(22)

This step can be omitted by using the *sigma points* from equation (18) to enhance the computational efficiency at the cost of performance [3]. The observation function is used to transform the *sigma points* into predicted measurements, vector.

(23)

The average of is the predicted measurement at time-step :

(24)

Equations (23 and 24) are used to estimate the covariance of predicted measurement

(25)

where takes into account the measurement noise.

Next, we estimate the cross covariance between and

(26)

Finally, the measurement at the time-step is used to update the state vector and its covariance

(27)

where

(28)

The *a posteriori* observation is recovered by .

We calculate the AIC measure for the two models used in Figure 6 using the following equations [35]

(29)

Where is the total number of parameters in the model, is the total number of data samples ( for examples in Figure 6), and is the residual sum of squares. The model that includes ion concentration dynamics has four extra parameters, , , , and . Therefore, we take =0 and 4 for the models without and with ion concentrations dynamics respectively.

All simulations were carried out using MATLAB on 24 multi-core Mac Pro computer. The MATLAB code for the results is archived at ModelDB (http://senselab.med.yale.edu/modeldb/default.asp) and can also be provided by the authors upon request.

Estimates of remaining variables for the INs shown in Figure 7. (A) intracellular *Ca ^{2+}* concentration (arbitrary units), (B)

(0.40 MB TIF)

Click here for additional data file.^{(393K, tif)}

Estimates of remaining variables for the PCs shown in Figure 8A, B. intracellular *Ca ^{2+}* concentration (arbitrary units) (A) and gating variables,

(0.30 MB TIF)

Click here for additional data file.^{(295K, tif)}

Estimates of synaptic variables for PCs and INs shown in Figure 8A–D. Synaptic variables, *s _{p}* (A),

(0.21 MB TIF)

Click here for additional data file.^{(202K, tif)}

We extend our heartfelt thanks to Jokubas Ziburkus for his constructive comments on the manuscript and generously providing us with access to the experimental data.

The authors have declared that no competing interests exist.

This study was supported by NIH Grants No. R01MH50006 and No. K02MH01493, the Pennsylvania Keystone Innovation Zone Program and Tobacco Settlement, and a grant from the National Academies - Keck Futures Initiative. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

1. Kalnay E. Atmospheric modeling, data assimilation, and predictability. 2003. Cambridge University Press, New York.

2. Kalman RE. A new approach to linear filtering and prediction problems. Trans ASME J Basic Eng. 1960;82:35–45.

3. Simon D. Optimal state estimation. Wiley-Interscience; 2006.

4. Baek SJ, Hunt BR, Kalnay E, Ott E, Szunyogh I. Local ensemble Kalman filtering in the presence of model bias. Tellus A. 2006;58:293–306.

5. Yang SC, Baker D, Li H, Cordes K, Huff M, et al. Data assimilation as synchronization of truth and model: Experiments with the three-variable Lorenz system. J Atmos Sci. 2006;63:2340–2354.

6. Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. J Geophys Res. 1994;99:10143–10162.

7. Evensen G, van Leeuwen PJ. An ensemble Kalman smoother for nonlinear dynamics. Mon Weather Rev. 2000;128:1852–1867.

8. Julier SJ, Uhlmann JK. A consistent, debiased method for converting between polar and cartesian coordinate systems. P SPIE. 1997;3068:110–121.

9. Julier SJ, Uhlmann JK. A new extension of the kalman filter to nonlinear systems. P SPIE. 1997;3068:182–193.

10. Cressman JR, Jr, Ullah G, Ziburkus J, Schiff SJ, Barreto E. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. J Comput Neurosci. 2009;26:159–170. [PMC free article] [PubMed]

11. Voss HU, Timmer J. Nonlinear dynamics system identification from uncertain and indirect measurements. Int J Bifurcat Chaos. 2004;14:1905–1933.

12. Schiff SJ, Sauer TD. Kalman filter control of a model of spatiotemporal cortical dynamics. J Neur Eng. 2008;5:1–8. [PMC free article] [PubMed]

13. Ullah G, Schiff SJ. Tracking and control of neuronal Hodgkin-Huxley dynamics. Phys Rev E. 2009;79:040901. [PMC free article] [PubMed]

14. Sauer TD, Schiff SJ. Data assimilation for heterogeneous networks: The consensus set. Phys Rev E. 2009;79:051909. [PMC free article] [PubMed]

15. Schiff SJ. Towards model-based control of parkinson's disease. Phil Trans Royal Soc A. 2010 in press. [PMC free article] [PubMed]

16. Schiff SJ, Cressman JR, Jr, Barreto E. Towards a dynamics of seizure mechanics. 2008. pp. 496–514. in Computational Neuroscience in Epilepsy (Academic Press), London:

17. Lehnertz K, Mormann F, Osterhage H, Muller A, Prusseit J, et al. State-of-the-art of seizure prediction. J Clin Neurophysiol. 2007;24:147–153. [PubMed]

18. Mormann F, Andrzejak RG, Elger CE, Lehnertz K. Seizure prediction: the long and winding road. Brain. 2007;130:314–333. [PubMed]

19. Mormann F, Kreuz T, Rieke C, Andrzejak RG, Kraskov A, et al. On the predictability of epileptic seizures. J Clin Neurophysiol. 2005;116:569–587. [PubMed]

20. Shu Y, Hasenstaub A, McCormick DA. Turning on and off recurrent balanced cortical activity. Nature. 2003;423:288–293. [PubMed]

21. Chub N, Mentis GZ, O'Donovan MJ. Chloride-sensitive MEQ fluorescence in chick embryo motoneurons following manipulations of chloride and during spontaneous network activity. J Neurophysiol. 2006;95:323–330. [PubMed]

22. Huang X, Troy WC, Yang Q, Ma H, Laing CR, et al. Spiral waves in disinhibited mammalian neocortex. J Neurosci. 2004;24:9897–9902. [PMC free article] [PubMed]

23. Ziburkus J, Cressman JR, Jr, Barreto E, Schiff SJ. Interneuron and pyramidal cell interplay during in vitro seizure-like events. J Neurophysiol. 2006;95:3948–3954. [PMC free article] [PubMed]

24. Schiff SJ, Sauer TD, Kumar R, Weinstein SL. Neuronal spatiotemporal pattern discrimination: the dynamical evolution of seizures. Neuroimage. 2005;28:1043–1055. [PMC free article] [PubMed]

25. Ullah G, Cressman JR, Jr, Barreto E, Schiff SJ. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: II. Network and glial dynamics. J Comput Neurosci. 2009;26:171–183. [PMC free article] [PubMed]

26. Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ. Slow state transitions of sustained neural oscillations by activity-dependent modulation of intrinsic excitability. J Neurophysiol. 2004;92:1116–1132. [PMC free article] [PubMed]

27. Frohlich F, Timofeev I, Sejnowski TJ, Bazhenov M. Extracellular potassium dynamics and epileptogenesis. 2008. pp. 419–439. In: Computational Neuroscience in Epilepsy (Academic Press), London:

28. Mitra P, Bokil H. Observed brain dynamics. 2007. Oxford University Press, USA.

29. Toth Z, Peña M. Data assimilation and numerical forecasting with imperfect models: The mapping paradigm. Physica D. 2007;230:146–158.

30. Ullah G, Schiff SJ. Models of epilepsy. Scholarpedia. 2009;4:1409.

31. Somjen GG. Ions in the brain: normal function, seizures, and stroke. 2004. Oxford University Press, USA.

32. Traynelis SF, Dingledine R. Potassium-induced spontaneous electrographic seizures in the rat hippocampal slice. J Neurophysiol. 1988;59:259–276. [PubMed]

33. Jensen MS, Yaari Y. Role of intrinsic burst firing, potassium accumulation, and electrical coupling in the elevated potassium model of hippocampal epilepsy. J Neurophysiol. 1997;77:1224–1233. [PubMed]

34. Bikson M, Hahn PJ, Fox JE, Jefferys JGR. Depolarization block of neurons during maintenance of electrographic seizures. J Neurophysiol. 2003;90:2402–2408. [PubMed]

35. Akaike H. A new look at the statistical identification model. IEEE T Automat Contr. 1974;19:716–723.

36. Trevelyan AJ, Sussillo D, Watson BO, Yuste R. Modular propagation of epileptiform activity: Evidence for an inhibitory veto in neocortex. J Neurosci. 2006;26:12447–12455. [PubMed]

37. Brown EN, Nguyen DP, Frank LM, Wilson MA, Solo V. An analysis of neural receptive field plasticity by point process adaptive filtering. P Natl Acad Sci USA. 2001;98:12261–12266. [PubMed]

38. Barbieri R, Wilson MA, Frank LM, Brown EN. An analysis of hippocampal spatio-temporal representations using a Bayesian algorithm for neural spike train decoding. IEEE T Neur Sys Reh. 2005;13:131–136. [PubMed]

39. Smith AC, Brown EN. Estimating a state-space model from point process observations. Neural Comp. 2003;15:965–991. [PubMed]

40. Srinivasan L, Brown EN. A state-space framework for movement control to dynamic goals through brain-driven interfaces. IEEE T Bio-med Eng. 2007;54:526–535. [PubMed]

41. Srinivasan L, Eden UT, Mitter SK, Brown EN. General-purpose filter design for neural prosthetic devices. J Neurophysiol. 2007;98:2456–2475. [PubMed]

42. Smith AC, Wirth S, Suzuki WA, Brown EN. Bayesian analysis of interleaved learning and response bias in behavioral experiments. J Neurophysiol. 2007;97:2516–2524. [PubMed]

43. Wu W, Gao Y, Bienenstock E, Donoghue JP, Black MJ. Bayesian population decoding of motor cortical activity using a Kalman filter. Neural Comp. 2006;18:80–118. [PubMed]

44. Ergun A, Barbieri R, Eden UT, Wilson MA, Brown EN. Construction of point process adaptive filter algorithms for neural systems using sequential Monte Carlo methods. IEEE T Bio-med Eng. 2007;54:419–428. [PubMed]

45. Li Z, O'Doherty JE, Hanson TL, Lebedev MA, Henriquez CS, et al. Unscented Kalman Filter for Brain-Machine Interfaces. PLoS One. 2009;4:e6243. [PMC free article] [PubMed]

46. Hunt BR, Kostelich EJ, Szunyogh I. Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter. Physica D. 2007;230:112–126.

47. Ott E, Hunt BR, Szunyogh I, Zimin AV, Kostelich EJ, et al. Estimating the state of large spatio-temporally chaotic systems. Phys Lett A. 2004;330:365–370.

48. Baek SJ, Hunt BR, Kalnay E, Ott E, Szunyogh I. Local ensemble Kalman filtering in the presence of model bias. Tellus A. 2006;58:293–306.

49. Spiller ET, Budhiraja A, Ide K, Jones CKRT. Modified particle filter methods for assimilating Lagrangian data into a point-vortex model. Physica D. 2008;237:1498–1506.

50. Salman H, Kuznetsov L, Jones C, Ide K. A method for assimilating lagrangian data into a shallow-water-equation ocean model. Mon Weather Rev. 2006;134:1081–1101.

51. Cornick M, Hunt BR, Ott E, Kurtuldu H, Schatz MF. State and parameter estimation of spatiotemporally chaotic systems illustrated by an application to Rayleigh–Bénard convection. Chaos. 2009;19:013108. [PubMed]

52. Paninski L, Ahmadian Y, Ferreira DG, Koyama S, Rahnama Rad K, et al. A new look at state-space models for neural data. J Comput Neurosci. 2009:1–20.

53. Huys QJM, Paninski L. Smoothing of, and parameter estimation from, noisy biophysical recordings. PLoS Comput Biol. 2009;5:e1000379. [PMC free article] [PubMed]

54. Deng B, Wang J, Che Y. A combined method to estimate parameters of neuron from a heavily noise-corrupted time series of active potential. Chaos. 2009;19:015105. [PubMed]

55. Charney JG. Dynamic forecasting by numerical process. 1951. pp. 470–482. in Compendium of Meteorology, American Meteorological Society, Boston:

56. Hodgkin AL, Huxely A. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117:500–544. [PubMed]

57. Lauger P. Electrogenic ion pumps. 1991. Distinguished lecture series of the Society of General Physiologists, Vol 5, Sinauer Associates Inc MA, USA.

58. Kager H, Wadman WJ, Somjen GG. Seizure-like afterdischarges simulated in a model neuron. J Comput Neurosci. 2007;22:105–128. [PubMed]

59. Paulson OB, Newman EA. Does the release of potassium from astrocyte endfeet regulate cerebral blood flow? Science. 1987;237:896–898. [PMC free article] [PubMed]

60. Kuschinsky W, Wahl M, Bosse O, Thurau K. The dependency of the pial arterial and arteriolar resistance on the perivascular h+ and k+ conconcentrations. a micropuncture study. Eur Neurol. 1972;6:92–95. [PubMed]

61. McCulloch J, Edvinsson L, Watt P. Comparison of the effects of potassium and ph on the calibre of cerebral veins and arteries. Pflugers Arch. 1982;393:95–98. [PubMed]

62. Fisher RS, Pedley TA, Prince DA. Kinetics of potassium movement in norman cortex. Brain Res. 1976;101:223–237. [PubMed]

63. Scharrer E. The blood vessels of the nervous tissue. Quart Rev Biol. 1944;19:308–318.

64. McBain CJ, Traynelis SF, Dingledine R. Regional variation of extracellular space in the hippocampus. Science. 1990;249:674–677. [PubMed]

65. Gutkin BS, Laing CR, Colby CL, Chow CC, Ermentrout GB. Turning on and off with excitation: the role of spike-timing asynchrony and synchrony in sustained neural activity. J Comput Neurosci. 2001;11:121–134. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

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. |