Bayesian Belief Updating of Spatiotemporal Seizure Dynamics

Bayesian Belief Updating of Spatiotemporal Dynamics Gerald Cooray 1 2 Richard Rosch 2 3 T orsten Baldeweg 3 Louis Lemieux 1 Karl Friston 2 Biswa Sengupta 4 5 6 Abstract Epileptic seizure activity shows complicated dy- namics in both space and time. T o understand the ev olution and propagation of seizures spa- tially extended sets of data need to be analysed. W e hav e previously described an efficient filter - ing scheme using variational Laplace that can be used in the Dynamic Causal Modelling (DCM) framew ork (Friston et al., 2003) to estimate the temporal dynamics of seizures recorded using ei- ther inv asi ve or non-in vasi ve electrical record- ings (EEG/ECoG). Spatiotemporal dynamics are modelled using a partial differential equation – in contrast to the ordinary differential equation used in our previous w ork on temporal estimation of seizure dynamics (Cooray et al., 2016). W e provide the requisite theoretical background for the method and test the ensuing scheme on simu- lated seizure acti vity data and empirical in vasi ve ECoG data. The method provides a framework to assimilate the spatial and temporal dynamics of seizure activity , an aspect of great physiologi- cal and clinical importance. 1. Introduction Epilepsy is a chronic disorder characterised by heteroge- neous and dynamic pathophysiological processes that lead to an altered balance between excitatory and inhibitory in- fluences at the cortical le vel. Electrocorticography (ECoG) recordings use a grid of electrodes to cov er a cortical area of clinical importance. This methodology allo ws for a sam- pling of epileptic seizure acti vity . Quantitati ve analysis of the temporal aspect of cortical activity using neural mass 1 Department of Clinical and Experimental Epilepsy , Institute of Neurology , University College London, UK 2 W ellcome T rust Centre for Neuroimaging, Institute of Neurology , Uni versity Col- lege London, UK 3 Dev elopmental Neuroscience Programme, In- stitute of Child Health, Uni versity College London, UK 4 Dept. of Bioengineering, Imperial College London 5 Dept. of Engineering, Univ ersity of Cambridge 6 Cortexica V ision Systems, UK. Corre- spondence to: Biswa Sengupta < b .sengupta@imperial.ac.uk > . ICML 2017 Time Series W orkshop, Sydney , Australia, 2017. Copyright 2017 by the author(s). (mean-field) models was introduced by W ilson and Co wan (Jansen & Rit, 1995; Wilson & Cow an, 1972; 1973). V ari- ation in the parameters of neural mass models hav e subse- quently been used to model seizure acti vity , and inferences about these changes can be made using Bayesian statistics; including stochastic filtering or genetic algorithms (Blenk- insop et al., 2012; Freestone et al., 2014; Nev ado-Holgado et al., 2012; Schiff & Sauer, 2008; Ullah & Schiff, 2010; W endling et al., 2005; Sengupta et al., 2015a;b). The in- herent assumptions of neural mass models preclude an in depth analysis of spatiotemporal dynamics, howe ver , other neuronal models, such as neural field models, have explic- itly included the spatial aspects of neuronal acti vity (Moran et al., 2013; Pinotsis et al., 2012; Schiff & Sauer, 2008). The comparati vely high computational comple xity of these models makes inference under such models very dif ficult. Sev eral studies have suggested that during seizures there are changes in intra- and extracellular factors such as electrolytes, metabolites and neurotransmitters. Mediated through the local interaction between glial and neuronal cells many of these factors will change the intrinsic dy- namics of the cortex (Fr ¨ ohlich et al., 2010; Ullah & Schiff, 2010; W ei et al., 2014). W e argue that a simplified spa- tiotemporal field can model the temporal and spatial be- haviour of seizure activity . T o contain model complex- ity and allow inference, we decouple fast temporal activ- ity from slower spatiotemporal fluctuations. This adia- batic assumption respects the biology based on the diffu- sion or transport of extracellular or intracellular factors as described above and their interaction with the synaptic con- nectivity of cortical neurons, which are assumed to produce fast temporal acti vity . The model is formulated within the Dynamic Causal Mod- elling (DCM) framework (Friston et al., 2003) – a model comparison and averaging frame work for Bayesian infer- ence of hierarchical yet dynamical generativ e models. This paper introduces a novel application of dynamic causal modelling based upon a hierarchical DCM, with first (fast) lev el spectral activity and second (slow) lev el spatial dy- namics. Bayesian Belief Updating 2. Methods T o construct an in vertible generative model, we partition our model into hidden states/parameters that are not explic- itly (conditionally) dependent on spatial location ( x, θ i , θ e ) and parameters that depend on spatial location ( θ sp ) . Fur- thermore, we assume that the temporal fluctuations of the spatially dependent subset will be several orders slower than the temporal dynamics of the hidden states/parameters that are conditionally independent of spatial location. P a- rameters not governed by spatial dynamics are sampled from a Gaussian distribution. For reasons of numerical efficienc y , we introduce the simplifying assumption that transport of electrolytes through the extracellular medium or through glial cells is via passiv e diffusion. At this point we appeal to the assumption that the dynam- ics of the neuronal activity ( x ) compared to the spatially varying parameters ( θ sp ) are at least sev eral orders of mag- nitude faster . This allows us to estimate neuronal activity as the steady state activity of a neural mass governing x . W e can then formulate the expected spectrum of the neu- ronal activity as a function of the parameters ( θ i , θ sp , θ e ) . This approach has been adopted in a similar setting (Cooray et al., 2016; Moran et al., 2011). Our model of the spectra measured during discrete time in- tervals from the local field potentials is given by the follow- ing measurement model and a partial differential equation (Eq. 1): y k = H ( θ i ( t k ) , θ sp ( t k ) , θ e ) + r k dθ sp dt = O 2 θ sp + q ( t ) (1) In the follo wing equations, sampling error r k ∼ N (0 , R k ) is sampled from a Gaussian distribution and q ( t ) is defined as a white noise process (Eq. 2). The diffusion equation is infinite dimensional, making it difficult to solve analyti- cally and numerically . Howe ver , we approximate the par- tial differential equation using a finite set of eigenfunctions, φ i ( x ) , and eigenv alues, λ i . This results in the following set of equations (S ¨ arkk ¨ a et al., 2012; Solin et al., 2012): y k = H ( θ i ( x, t k ) , θ sp ( x, t k ) , θ e ( x, t k )) + r k θ sp ( t ) ≈ X c i ( t k ) φ i ( x ) e − λ i ( t − t k ) c i ( t k ) = c i ( t k − 1 ) e − λ i ( t k − t k − 1 ) + q i ( t k − 1 ) (2) In this equation q i ( t k ) ∼ N (0 , Q i,k ) is sampled from a Gaussian distribution. A similar idea has previously been presented in DCM of distributed electromagnetic responses where the cortical activity was parameterised using a set Figure 1. Left pane: Simulated fluctuation in parameter . Right pane: Results of in version of simulated data (not shown). of local standing-wa ves of neural field models (Daunizeau et al., 2009). The data generated by this model comprised of windowed spectral data y k,ij from the i th row and the j th column of any electrode array used to measure the electrographic ac- tivity sampled at time t k . W e estimate the parameters of our model in each time window using v ariational Laplace (MacKay, 2002). The values in the next time window are predicted using Eq. 2, given the previous values. W e can then estimate an approximate solution to the above stochas- tic model using Bayesian belief updating as has been previ- ously described (Cooray et al., 2016). For technical details please refer to Section A-D. 3. Results W e simulated ECoG data with two-dimensional spatial dy- namics: giving an excellent fit to the spectral activity (ex- plained variation > 0.99; Figure 1). Furthermore, the in- ferred fluctuations of the excitatory parameter were nearly identical to the simulated dynamics. W e also estimated a spatially dependent excitatory gain parameter for epileptic seizure activity using recordings from two patients with focal cortical dysplasia and refrac- tory epilepsy (Figure 2). Subdural recordings were made from a 32-contact and a 20-contact platinum array (Ad- T ech) placed over the lesion. W e used two seizures free of artefacts from each data set for modelling the averaged in- duced spectral activity . The model also estimated parame- ters (without spatiotemporal dynamics) of intrinsic connec- tivity within the cortical column that were sampled from a white noise process. W e allo wed this random variation dur- ing Bayesian belief updating to accommodate unmodelled changes during the seizure activity . There was a good fit between the data and the predicted activity (explained vari- ation > 0.85). Bayesian Belief Updating Figure 2. Estimated cortical excitability at dif ferent time points Spatial dynamics controlling cortical excitability were in- ferred from the data, enabling parameters (excitatory gain) to be estimated at each point within the cortical grid ov er time. The parameter increased to a maximum tow ards the end of the seizure (a) and then quickly dissipates during termination of the seizure (b) (Figure 2). 4. Discussion W e present a framework for analysing spatiotemporal dy- namics, in our specific case such a dynamics emerges from epileptic seizure activity . W e extended an inference scheme previously described for inference of the temporal structure of seizure acti vity (Cooray et al., 2016; Cooray et al., 2017) – we achiev ed this objectiv e using a field model to incor- porate spatial dynamics. Including spatial processes allo w us to model the epileptic spread of seizure activity , which is one of the key physiological processes underlying a fo- cal seizure. A deeper understanding of the spatial spread of focal seizures ov er the cortex could be used for curtailing the spread of a seizure; limiting its effect on the patient by either surgical or other procedures such as direct or induced electrical currents delivered by e.g. transcranial magnetic stimulation or deep brain stimulation (Engel, 1993; Roten- berg et al., 2008). The current frame work for the analysis of seizure activity assumes that the spatial dynamics can be estimated using a finite set of eigenfunctions of the underlying partial dif- ferential equation (modelling the spatial dynamics). This assumption is only v alid for relati vely simple spatiotempo- ral dynamics (partial differential equations with analytical solutions) on simple representations of the cortical surface (such as a plane or spherical surface) (Ev ans, 2010). If any of these assumptions are violated, as in the case of realistic geometries obtained using MRI or more realistic descrip- tions of the spatiotemporal dynamics, the forward model would require re-formulation using numerical approxima- tions. One way forward is by reformulating our scheme under a finite element framework using the approach de- scribed in Sengupta & Friston (2017); although this re- quires e xtensiv e computational expenditure. The real v alue of this work, therefore, lies in the potential to fit the pa- rameters of spatiotemporal dynamics of generativ e models to individual patients. This should allow us to understand individual differences that would not only guide epilepto- genic resections but also be used for prediction of epilepsy surgery outcomes. In conclusion, we present a straightforward w ay of extend- ing an established Bayesian belief update/filtering scheme to include the spatial spread of seizure activity . This re- quires several (plausible) assumptions regarding the be- haviour of cortical columns, the geometry of the cortex and a separation of time scales of f ast neuronal activity and slower fluctuations in cortical gain. W e suggest that these assumptions have some experimental evidence but that the method can be modified if these assumptions are shown to be inv alid or too inaccurate in the future. In particular , by comparing spatiotemporal models of different complexity (e.g. dif fusion models versus complex reaction-diffusion models with threshold dynamics), using their model evi- dence, the above framework could be used to address im- portant questions about how seizure activity spreads in pa- tients. References Blenkinsop, Alex, V alentin, Antonio, Richardson, Mark P , and T erry , John R. The dynamic ev olution of focal-onset epilepsies–combining theoretical and clinical observ ations. Eur opean Journal of Neur oscience , 36(2):2188–2200, 2012. Cooray , Gerald K, Sengupta, Biswa, Douglas, Pamela K, and Friston, Karl. Dynamic causal modelling of electrographic seizure activity using Bayesian belief updating. Neur oImage , 125:1142–1154, 2016. Cooray, G.K, Rosch, R., Balde weg, T ., Lemieux, L., Friston, K., and Sengupta, B. Bayesian Belief Updating of Spatiotemporal Seizure Dynamics. ArXiv e-prints , 2017. Daunizeau, Jean, Kiebel, Stefan J, and Friston, Karl J. Dynamic causal modelling of distributed electromagnetic responses. Neur oImage , 47(2):590–601, 2009. Engel, J. Surgical tr eatment of the epilepsies . Raven Press, 1993. Evans, CE. P artial Differ ential Equations . American Mathemat- ical Society , 2010. Freestone, Dean R, Karoly , Philippa J, Ne ˇ si ´ c, Dragan, Aram, Parham, Cook, Mark J, and Grayden, Da vid B. Estimation of effecti ve connecti vity via data-driven neural modeling. F r on- tiers in neur oscience , 8, 2014. Friston, Karl J, Harrison, Lee, and Penny , Will. Dynamic causal modelling. Neuroimag e , 19(4):1273–1302, 2003. Fr ¨ ohlich, Flavio, Sejnowski, T errence J, and Bazhenov , Maxim. Network bistability mediates spontaneous transitions between normal and pathological brain states. Journal of Neur oscience , 30(32):10734–10743, 2010. Bayesian Belief Updating Jansen, Ben H and Rit, V incent G. Electroencephalogram and visual ev oked potential generation in a mathematical model of coupled cortical columns. Biological cybernetics , 73(4):357– 366, 1995. MacKay , Da vid J. C. Information Theory , Infer ence & Learning Algorithms . Cambridge Univ ersity Press, New Y ork, NY , USA, 2002. Moran, Rosalyn, Pinotsis, Dimitris A, and Friston, Karl. Neural masses and fields in dynamic causal modeling. F r ontiers in computational neur oscience , 7, 2013. Moran, Rosalyn J, Stephan, Klaas E, Dolan, Raymond J, and Fris- ton, Karl J. Consistent spectral predictors for dynamic causal models of steady-state responses. Neuroima ge , 55(4):1694– 1708, 2011. Nev ado-Holgado, Alejo J, Marten, Frank, Richardson, Mark P , and T erry , John R. Characterising the dynamics of EEG wa ve- forms as the path through parameter space of a neural mass model: application to epilepsy seizure e volution. Neuroima ge , 59(3):2374–2392, 2012. Pinotsis, Dimitris A, Moran, Rosalyn J, and Friston, Karl J. Dy- namic causal modeling with neural fields. Neur oimage , 59(2): 1261–1274, 2012. Rotenberg, Ale xander , Muller, Paul, Birnbaum, Daniel, Harring- ton, Michael, Riviello, James J, Pascual-Leone, Alvaro, and Jensen, Frances E. Seizure suppression by EEG-guided repeti- tiv e transcranial magnetic stimulation in the rat. Clinical Neu- r ophysiology , 119(12):2697–2702, 2008. S ¨ arkk ¨ a, Simo, Solin, Arno, Nummenmaa, Aapo, V ehtari, Aki, Auranen, T oni, V anni, Simo, and Lin, Fa-Hsuan. Dynamic retrospectiv e filtering of physiological noise in BOLD fMRI: DRIFTER. NeuroImag e , 60(2):1517–1527, 2012. Schiff, Ste ven J and Sauer , T im. Kalman filter control of a model of spatiotemporal cortical dynamics. BMC Neuroscience , 9(1): O1, 2008. Sengupta, B., Friston, K. J., and Penny , W . D. Gradient-based MCMC samplers for dynamic causal modelling. Neuroima ge , 2015a. Sengupta, B., Friston, K. J., and Penny , W . D. Gradient-free MCMC methods for dynamic causal modelling. Neuroimag e , 112:375–81, 2015b. Sengupta, Biswa and Friston, Karl. Sentient self-organization: Minimal dynamics and circular causality . arXiv preprint arXiv:1705.08265 , 2017. Solin, Arno et al. Hilbert space methods in infinite-dimensional Kalman filtering. 2012. Ullah, Ghanim and Schif f, Stev en J. Assimilating seizure dynam- ics. PLoS computational biology , 6(5):e1000776, 2010. W ei, Y ina, Ullah, Ghanim, Ingram, Justin, and Schiff, Stev en J. Oxygen and seizure dynamics: II. Computational modeling. Journal of neur ophysiolo gy , 112(2):213–223, 2014. W endling, Fabrice, Hernandez, Alfredo, Bellanger , Jean-Jacques, Chauvel, Patrick, and Bartolomei, Fabrice. Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intrac erebral EEG. J ournal of Clinical Neur ophysiology , 22(5):343, 2005. W ilson, Hugh R and Cowan, Jack D. Excitatory and inhibitory interactions in localized populations of model neurons. Bio- physical journal , 12(1):1–24, 1972. W ilson, Hugh R and Cowan, Jack D. A mathematical theory of the functional dynamics of cortical and thalamic nervous tis- sue. Biological Cybernetics , 13(2):55–80, 1973. Bayesian Belief Updating A ppendices A. Canonical mean-field equations The cortical columns described in this note were modelled using a neural mass (mean-field) model generating electro- graphic seizure activity . The canonical cortical microcircuit (CMC) is comprised of four subpopulations of neurons cor- responding to superficial and deep pyramidal, excitatory , and inhibitory cells (Moran et al., 2013). These cell populations are connected using ten inhibitory and e xcitatory connections. Afferent connections dri ve the excitatory granular cells and efferent connections derive from the superficial pyramidal cells. The equations gov erning the four mean-fields are giv en by , ¨ x e + 2 ˙ x e T e + x e T 2 e = − g 1 s ( x e ) − g 3 s ( x i ) − g 2 s ( x sp ) ¨ x i + 2 ˙ x i T i + x i T 2 i = θ sp g 5 s ( x e ) + θ sp g 6 s ( x dp ) − g 4 s ( x i ) ¨ x sp + 2 ˙ x sp T sp + x sp T 2 sp = θ sp g 8 s ( x e ) − g 7 s ( x sp ) ¨ x dp + 2 ˙ x dp T dp + x dp T 2 dp = − g 10 s ( x dp ) − g 9 s ( x i ) (3) θ sp is the parameter with spatiotemporal dynamics. Parameter g 7 was sampled from a Gaussian distribution during Bayesian filtering. All other parameters were kept constant during the in version. B. Solution of parameter field equations The spatio-temporal (on a single spatial dimension) parameter is modelled as a heat equation with time dependent boundary conditions, u ( x, t ) : parameter with spatiotemporal variation α : diffusion coef ficient φ 1 / 0 ( t ) : boundary condition (time dependent) f ( x ) : initial condition L : length of one dimensional domain ov er which u varies ω ( x ) : time independent particular solution v ( x, t ) : homogenous solution to diffusion equation S ( x, t ) : auxillary variable u t = α 2 u xx u (0 , t ) = φ 0 , u ( L, t ) = φ 1 u ( x, 0) = f ( x ) ω = φ 0 + x L ( φ 1 − φ 0 ) u = ω + v u t = ω t + v t = α 2 v xx v t = α 2 v xx − ω t v (0 , t ) = v ( L, t ) = 0 Bayesian Belief Updating v ( x, 0) = f ( x ) − ω ( x, 0) S ( x, t ) = − ˙ ω (4) The solution to the diffusion equation is separated into time and space dependent functions; i.e., by using separation of variables. The full solution can then be written as an infinite series of eigenfunctions (space dependent functions), to- gether with their time-dependent variation. The eigenfunctions of the diffusion equation satisfy an equiv alent ordinary differential equation. For the one-dimensional diffusion equation with Dirichlet boundary conditions (zero along bound- aries), the eigenfunctions can be written as sine functions. For the two-dimensional case with rotational symmetry , a linear combination of first and second order Bessel functions are used, λ n : eigen v alue of ordinary differential equation S n ( t ) : coefficient of the eigenfunction sinλ n x in the eigenfunction e xpansion of the auxiliary variable S ( x, t ) v n ( t ) : coefficient of the eigenfunction sinλ n x in the eigenfunction e xpansion of v ( x, t ) c n : coef ficient of the eigenfunction sinλ n x in the eigenfunction e xpansion of the initial condition of v ( x, t ) , i.e., v ( x, 0) S ( x, t ) = ∞ X n =1 S n ( t ) sinλ n x S n ( t ) = − 2 L Z L 0 dω dt sinλ n xdx λ n = nπ L v ( x, t ) = ∞ X n =1 v n ( t ) sinλ n x v ( x, 0) = ∞ X n =1 c n sinλ n x c n = 2 L Z L 0 ( f ( x ) − ω ( x, 0)) sinλ n xdx u ( x, t ) = ω ( x, t ) + ∞ X n =1  Z t 0 S n ( t ) exp  α 2 λ 2 n τ  dτ + c n  exp  − α 2 λ 2 n t  sinλ n x (5) C. Bayesian belief updating The following equations operationalize the Bayesian belief updating of the parameters as data is inv erted sequentially across the windowed data, y i , see (Cooray et al., 2016; Cooray et al., 2017) for more details. The priors on the n eigen- function coefficients of parameter θ sp in each window is gi ven by , µ i : posterior mean v ector of θ sp in the ith window i.e., ( c 1 , c 2 , . . . , c n ) i λ : eigen values, i.e. ( λ 1 , λ 2 , . . . , λ n ) ∆ : time interval between windo ws Q i : posterior co variance matrix of θ sp in the ith window R : volatility cov ariance matrix of θ sp p ( θ sp | y i , . . . , y 1 ) = prior probability density of θ sp in the i+1th window p ( θ sp | y i , . . . , y 1 ) = N  µ i e − λ ∆ , Q i + R  (6) Bayesian Belief Updating Using a variational Bayes formulation, we update this prior probability distrib ution using the data in window i + 1 . D. P arameterizations of field models Neural field u ( x, t ) : neural field representing acti vity of a population of neurons at location x at time t w ( y ) : the strength of connections between neurons separated by a distance y f ( u ) : firing rate function, u t ( x, t ) = − u + Z ∞ −∞ w ( y ) f ( u ( x − y , t )) dy (7) Parameter field Note that the spatial variation of the neural field is due to the time independent variable w . As this term is con volv ed ov er space with the firing rate it becomes difficult to separate the spatial and the temporal dynamics. In the case of the parameter field it is possible to have a slowly diffusing parameter (diffusion rate dependent on α ) that interacts as though it is constant with respect to the dynamics of the faster neural mean fields (rate dependent on k ). u ( x, t ) : neural field representing activity of a population of neurons at location x at time t f ( u ) : firing rate function θ ( x, t ) : parameter field k : rate coefficient of neural mass α : diffusion coef ficient u tt ( x, t ) − 2 k u t ( x, t ) − k 2 u ( x, t ) = f ( u ( x, t ) , θ ( x, t )) θ t = α 2 θ xx (8)