Assimilating Seizure 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 ., 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 Akaikes 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 ., 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 ., 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 .
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 .
neuroscience/theoretical neuroscience, computational biology/computational neuroscience, neurological disorders/epilepsy