A data-driven framework for neural field modeling

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/49844161 A data-driven framework for neural field modeling Article in NeuroImage · February 2011 DOI: 10.1016/j.neuroimage.2011.02.027 · Source: PubMed CITATIONS READS 31 148 6 authors, including: Dean Robert Freestone Parham Aram 47 PUBLICATIONS 219 CITATIONS 17 PUBLICATIONS 64 CITATIONS University of Melbourne SEE PROFILE The University of Sheffield SEE PROFILE K. Scerri Visakan Kadirkamanathan 8 PUBLICATIONS 80 CITATIONS 198 PUBLICATIONS 2,500 CITATIONS University of Malta SEE PROFILE The University of Sheffield SEE PROFILE Some of the authors of this publication are also working on these related projects: Neurovista Study View project DESIGN OF A SMART INFLATED TORUS AND ANTENNA MEMBRANE View project All content following this page was uploaded by Parham Aram on 23 June 2015. The user has requested enhancement of the downloaded file. A Data-Driven Framework for Neural Field Modeling D. R. Freestone a,b,c,1,∗ , P. Aramd,1 , M. Dewarc,e , K. Scerrif , D. B. Graydena,b , V. Kadirkamanathand a Department of Electrical and Electronic Engineering, University of Melbourne, Melbourne, VIC, Australia b The Bionic Ear Institute, East Melbourne, VIC, Australia c Institute for Adaptive and Neural Computation, University of Edinburgh, Edinburgh, UK d Department of Automatic Control and Systems Engineering, University of Sheffield, Sheffield, UK e Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY, USA f Department of Systems and Control Engineering, University of Malta, Msida, MSD, Malta Abstract This paper presents a framework for creating neural field models from electrophysiological data. The Wilson and Cowan or Amari style neural field equations are used to form a parametric model, where the parameters are estimated from data. To illustrate the estimation framework, data is generated using the neural field equations incorporating modeled sensors enabling a comparison between the estimated and true parameters. To facilitate state and parameter estimation, we introduce a method to reduce the continuum neural field model using a basis function decomposition to form a finite-dimensional state-space model. Spatial frequency analysis methods are introduced that systematically specify the basis function configuration required to capture the dominant characteristics of the neural field. The estimation procedure consists of a two-stage iterative algorithm incorporating the unscented Rauch-Tung-Striebel smoother for state estimation and a least squares algorithm for parameter estimation. The results show that it is theoretically possible to reconstruct the neural field and estimate intracortical connectivity structure and synaptic dynamics with the proposed framework. Keywords: Neural Field Model, Nonlinear Estimation, Intracortical Connectivity, Nonlinear Dynamics 1. Introduction Generating physiologically plausible neural field models is of great importance for studying brain dynamics at the mesoscopic and macroscopic scales. While our understanding of the function of neurons is well developed, the overall behaviour of the brain’s mesoscopic and macroscopic scale dynamics remains largely theoretical. Understanding the brain at this level is extremely important since it is at this scale that pathologies such as epilepsy, Parkinson’s disease and schizophrenia are manifested. ∗ Corresponding author Email addresses: [email protected] (D. R. Freestone ), [email protected] (P. Aram), [email protected] (M. Dewar), [email protected] (K. Scerri), [email protected] (D. B. Grayden), [email protected] (V. Kadirkamanathan) 1 Authors contributed equally to this work. Preprint submitted to Neuroimage June 23, 2015 Mathematical neural field models provide insights into the underlying physics and dynamics of electroencephalography (EEG) and magnetoencephalography (MEG) (see Deco et al. (2008); David and Friston (2003) for recent reviews). These models have demonstrated possible mechanisms for the genesis of neural rhythms (such as the alpha and gamma rhythms) (Liley et al., 1999; Rennie et al., 2000), epileptic seizure generation (Lopes Da Silva et al., 2003; Suffczynski et al., 2004; Wendling et al., 2005) and insights into other pathologies (Moran et al., 2008; Schiff, 2009) that would be difficult to gain from experimental data alone. Unfortunately, the use of these models in the clinic has been limited, since they are constructed for “general” brain dynamics whereas pathologies almost always have unique underlying patient-specific causes. Patient-specific data from electrophysiological recordings is readily available in the clinical setting, particularly from epilepsy surgery patients, suggesting an opportunity to make the patient-specific link to models of cortical dynamics. Furthermore, recent technological advances have driven an increased level of sophistication in recording techniques, with dramatic increases in spatial and temporal sampling (Brinkmann et al., 2009). However, the mesoscopic and macroscopic neural dynamic state is not directly observable in neurophysiological data, making predictions of the underlying physiology inherently difficult. For models to be clinically viable, they must be patient-specific. A possible approach to achieve this would be to fit a general continuum neural field model, like the Wilson and Cowan (1973) (WC) or Amari (1977) models, or a neural mass model like the Jansen and Rit (1995) model, to patient-specific data. Fitting the neural models to individuals is a highly non-trivial task. Recently, however, this task has been approached from a number of standpoints. The first paper on patient-specific modeling (to the authors’ best knowledge) came from Valdes et al. (1999), where they fit the neural mass model of Lopes Da Silva et al. (1976) and Zetterberg et al. (1978) to EEG data using the local linearization filter (Ozaki, 1993, 1994; Ozaki et al., 2000). This paper demonstrated for the first time the feasibility of assimilating data with neural models of EEG. Perhaps the most accepted estimation framework for neural mass models is dynamical causal modeling (DCM) (David and Friston, 2003; David et al., 2006), which has been proposed for studying evoked potential dynamics. This framework can be viewed as an extension to the work of Valdes et al. (1999) allowing for coupling of neural masses. Via a Bayesian inference scheme, DCM estimates the long range (cortico-cortical) connectivity structure between the specific isolated brain regions that best explains a given data set using the model of Jansen and Rit (1995). Data-driven neural modeling was extended to the continuum approximation by Galka et al. (2008), where they proposed an estimation framework based on a linear damped wave equation. Using a Kalman filter maximum likelihood framework, they improved on the low resolution electromagnetic tomography (LORETA) method for solving the inverse EEG problem. Another continuum neural field model-based estimator was developed by Daunizeau et al. (2009), as an extension to the DCM framework. They argue 2 that using a spatially extended continuum field model is superior for explaining cortical function than the neural mass DCM, since field models can explain a richer repertoire of dynamics such as traveling waves and bump solutions. Another recent approach estimates the parameters of a modified WC neural field model using an unscented Kalman filter (Schiff and Sauer, 2008). This work takes a systems theoretic approach to the neural estimation problem, successfully demonstrating that it is possible to perform state estimation and control on spatiotemporal neural fields. This marks the first step in what has the potential to revolutionize the treatment of many neurological diseases where therapeutic electrical stimulation is viable. For other valuable contributions to data-driven neural modeling, see Nunez (2000), Jirsa et al. (2002), and Robinson et al. (2004). We present an extension to the work of Schiff and Sauer (2008) by establishing a framework for estimating the state of the WC equations for larger scale (more space) systems via a systematic model reduction procedure. In addition, a new method is presented for estimating the connectivity structure and the synaptic dynamics. Until now, model-based estimation of local intracortical connectivity has not been reported in the literature (to the best of the authors’ knowledge). Our study also extends recent work which shows that it is possible to estimate local coupling of spatiotemporal systems using techniques from control systems theory and machine learning (Dewar et al., 2009). The key development of this previous work was to represent a spatiotemporal system as a standard state-space model, with the number of states independent of the number of observations (recording electrodes in this case). In addition, the appropriate model selection tools have been developed (Scerri et al., 2009) allowing for the application of the technique to neural fields. This paper extends the linear framework of Dewar et al. (2009) and Scerri et al. (2009) to the nonlinear case required for the neural field equations. Modeling the neural dynamics within this framework has a distinct advantage over the more standard multivariate auto-regressive (MVAR) models: the number of parameters to define the spatial connectivity is considerably smaller than the number of AR coefficients typically required to achieve the model complexity. In this paper, we demonstrate for the first time how an intracortical connectivity kernel can be inferred from data, based on a variant of the Wilson and Cowan (1973) neural field model. This work provides a fundamental link between the theoretical advances in neural field modeling and high resolution intracranial electrophysiological data. To illustrate the estimation framework, data is generated using the neural field equations incorporating modeled sensors enabling a comparison to be made between the estimated and true parameters. The paper proceeds by first describing the continuum neural field equations that are used as the cortical model. Then a finite-dimensional neural field model is derived. The model is reduced by approximating the neural field using a set of continuous basis functions, weighted by a finite dimensional state vector. The next section establishes conditions, using spatial frequency analysis, for both sensor and basis function spacing and width, such that the dominant dynamics of the neural field can be represented 3 by the reduced model. The state and parameter estimation procedure is described in the following section. The results for the spatial frequency analysis and parameter estimation are then presented. Finally, the implications and limitations of this framework are discussed along with planned future developments. 2. Methods 2.1. Neural Field Model Neural field models relate mean firing rates of pre-synaptic neural populations to mean post-synaptic membrane potentials. They are popular as they are parsimonious yet have a strong link with the underlying physiology. Each neural population represents a functional cortical processing unit, such as a column. The columnar organization of the cortex is continuous, where pyramidal cells are members of many columns. In general, cortical structure can be modeled in a physiologically plausible manner as being locally homogeneous (in short range intracortical connectivity) and heterogeneous (in long range cortico-cortical and corticothalamic connectivity) (Jirsa, 2009; Qubbaj and Jirsa, 2007). In certain regions of the cortex, each column is thought to be connected locally via symmetric short range local excitation with surround inhibition (Braitenberg and Schüz, 1998). For example, this structural organization is most studied in the visual system, where the surrounding inhibition effectively tunes a cortical column to a particular receptive visual field (Sullivan and De Sa, 2006). Neural field models are descriptive of a range of neurodynamics of the cortex such as evoked potentials, visual hallucinations and epileptic behaviour (David and Friston, 2003; Bressloff et al., 2001; Breakspear et al., 2006). Field models are also capable of generating complex spatial patterns of activity such as Turing patterns, spirals and traveling oscillations (Amari, 1977; Coombes, 2005; Coombes et al., 2007). It is an implicit assumption that the neural field model (in equation 12) provides an apt description of the cortical dynamics recorded from a specific subject. Although models of this form are capable of describing a variety of cortical dynamics, there will be without doubt a mismatch between the cortex and the model. Nevertheless, there is a sufficient volume of interesting results from theoretical studies using the WC field equations that warrants the data assimilation framework. An advantage in using a lumped-parameter field model over a more detailed mathematical description is that the myriad of parameters that might influence excitability, such as specific ion concentrations, are considered to be lumped into parameters (connectivity kernel coefficients) that effectively describe the net system gain. Therefore, there are less parameters to estimate. A major challenge in model-based data analysis for the mass action of the brain is to make models sufficiently detailed, such that the parameters are meaningful, but simple enough to yield clear insights that can be related to theoretical studies (from both neuroscience and engineering perspectives). 4 2.2. Integro-Difference Equation Neural Field Model The combination of modeling techniques in this paper leads to a large amount of notation, so a reference of the symbols used is provided in Table 1. The model relates the average number of action potentials g(r, t) Symbol Quantity Units Domain and indices Ω r t yt v(r, t) g(r, t) f (r, t) e(r, t) εt m(rn , r′ ) w(r, r′ ) h(t) ζ, ξ η(t) δ(t) D φ(r) xt ψ(r − r′ ) Ψ(r′ ) Γ Q() et C ν, ν c σν2 X W x̂ft − , x̂ft b x̂b− t , x̂t f− Pt , Ptf Ptb− , Ptb Mt Kt St q() Ŵ spatial domain spatial location time n.a. [mm, mm] s Model observation mean membrane potential weighted firing rate function firing rate function field disturbance, with covariance function γ(r − r′ ) observation noise, with covariance matrix Σε observation kernel where n is sensor index n = 1, .., ny (width σm ) spatial connectivity kernel post-synaptic response kernel inverse synaptic time constant & time constant parameter Heaviside step function Dirac-delta function temporal differential operator Reduced Model vector of nx Gaussian basis functions state vector at time t vector of nθ connectivity kernel basis functions decomposition connectivity matrix inner product of field basis functions state transition function state disturbance vector, with covariance Σe observation matrix Frequency Analysis spatial frequency and spatial cutoff frequency variance of Fourier transformed Gaussian basis function Estimation matrix of sigma vectors matrix of sigma vector weights (scaling parameter λ) forward prior and posterior state estimates backward prior and posterior state estimates forward prior and posterior covariance matrices backward prior and posterior covariance matrices cross-covariance matrix filter gain smoother gain nonlinear state function (used in LS estimator) vector of parameter estimates (used in LS estimator) mV mV spikes.s−1 spikes.s−1 mV mV n.a. n.a. mV s−1 , n.a. n.a. n.a. n.a. n.a. mV n.a. n.a. n.a. n.a. mV n.a. cycles/mm (cycles/mm)2 mV n.a. mV mV mV2 mV2 mV2 n.a. n.a. n.a. n.a. Table 1: Notation. The first column of the table gives the symbols for all the notation used in the paper. The second column provides a brief description of the quantity, and the third column provides the SI units where relevant. 5 arriving at time t and position r to the local post-synaptic membrane voltage v(r, t). The post-synaptic potentials generated at a neuronal population at location r by action potentials arriving from all other connected populations is described by v (r, t) = Z t h (t − t′ ) g (r, t′ ) dt′ . (1) −∞ The post-synaptic response kernel h(t) is described by h(t) = η(t) exp (−ζt), (2) where ζ = τ −1 , τ is the synaptic time constant and η(t) is the Heaviside step function. Non-local interactions between cortical populations at positions r and r′ are described by g (r, t) = Z w (r, r′ ) f (v (r′ , t)) dr′ , (3) Ω where f (·) is the firing rate function, w(·) is the spatial connectivity kernel and Ω is the spatial domain representing a cortical sheet or surface. The connectivity kernel is typically a “Mexican hat” or “wizard hat” function, which describes local excitation and surround inhibition. An example of the Mexican hat Connection Strength connectivity kernel, that is a weighted sum of Gaussians, used in this paper is shown in Fig. 1. To demon100 50 0 −50 0 −10 10 Space Figure 1: Mexican-hat connectivity kernel. The kernel is a composite of three Gaussian basis functions (dashed lines), representing short-range excitation, mid-range inhibition and long-range excitation. The kernel is rotationally symmetric (isotropic) about zero, hence a cross-section captures the important aspect of the connectivity kernel’s shape. The kernel decays asymptotically to zero. strate flexibility in the estimation algorithm and to model larger cortical regions, a third component of the kernel is introduced describing weak longer-range excitation. The exact shape of this kernel is assumed to vary across subjects, and hence needs to be inferred from subject-specific data. The firing rate of the presynaptic neurons is related to the post-synaptic membrane potential by the 6 sigmoidal activation function f (v (r′ , t)) = 1 . 1 + exp (ς (v0 − v (r′ , t))) (4) The parameter v0 describes the firing threshold of the neural populations and ς governs the slope of the sigmoid. By substituting equation 3 into equation 1 we get the spatiotemporal model v (r, t) = Z t ′ h (t − t ) −∞ Z w (r, r′ ) f (v (r′ , t′ )) dr′ dt′ . (5) Ω To obtain the standard integro-differential equation form of the model, we use the fact that the synaptic response kernel is a Green’s function of a linear differential equation defined by the differential operator D = d/dt + ζ. A Green’s function satisfies Dh (t) = δ (t) , (6) where δ(t) is the Dirac-delta function. Applying the differential operator D to equation 1 gives Dv (r, t) = D (h ∗ g) (r, t) (7) = (Dh ∗ g) (r, t) (8) = (δ (t) ∗ g) (r, t) (9) = g (r, t) (10) where ∗ denotes the convolution operator. This gives the standard form of the model dv (r, t) + ζv (r, t) = dt Z w (r, r′ ) f (v (r′ , t)) dr′ . (11) Ω To arrive at the integro-difference equation (IDE) form of the model, we discretize time using a first-order Euler method (see Appendix A for derivation) giving vt+Ts (r) = ξvt (r) + Ts Z w (r, r′ ) f (vt (r′ )) dr′ + et (r) , (12) Ω where Ts is the time step, ξ = 1−Ts ζ and et (r) is an i.i.d. disturbance such that et (r) ∼ GP(0, γ(r−r′ )). Here GP(0, γ(r − r′ )) denotes a zero mean Gaussian process with spatial covariance function γ(r − r′ ) (Rasmussen and Williams (2005)). The disturbance is added to account for model uncertainty and unmodeled inputs. To simplify the notation, the index of the future time sample, t + Ts , shall be referred to as t + 1 throughout the rest of the paper. The mapping between the membrane voltage and the electrophysiological data, denoted by yt , is modeled 7 using the observation function that incorporates sensors with a spatial extent by yt (rn ) = Z m (rn − r′ ) vt (r′ ) dr′ + εt (rn ), (13) Ω where m (rn − r′ ) is the observation kernel, rn defines the location of the electrodes in the field, n = 0, ..., ny − 1 indexes the sensors and εt (rn ) ∼ N (0, Σε ) denotes a multivariate normal distribution with mean zero and the covariance matrix Σε = σε2 I, where I is the identity matrix. Since we are considering intracranial measurements recorded directly from the surface of the cortex or within the brain, the lead field is not modeled by the observation equation. When considering the measurement of a point source in a conductive homogeneous field measured by a point sensor, the observation kernel will theoretically follow a Laplacian like shape, with an exponential decay in the amplitude of a measurement with an increasing distance between the source and sensor (Jackson, 1999). However, when considering a source or sensor with a spatial extent, the peak of the observation kernel becomes smoothed and is more similar to a Gaussian shape (Jackson, 1999). An experiment supporting this was performed and the results are shown in Appendix B. Therefore, in this study the observation kernel, m(r − r′ ), that governs the sensor pick-up geometry is defined by the Gaussian  (r − r′ )⊤ (r − r′ ) m (r − r ) = exp − , 2 σm  ′ (14) where σm sets the kernel width. The superscript ⊤ denotes the transpose operator. 2.3. Derivation of Finite Dimensional State-Space Model In order to implement standard estimation techniques, we use a decomposition of the field using a set of Gaussian basis functions. Decomposition allows a continuous field to be represented by a finitedimensional state vector. This facilitates application of standard nonlinear state estimation methods such as the unscented Kalman filter. The field decomposition is described by vt (r) ≈ φ⊤ (r) xt , (15) where xt is the state vector that scales the field basis functions φ(r). The field basis functions are described by ! (r − r′ )⊤ (r − r′ ) , φ (r − r ) = exp − σφ2 ′ (16) where σφ is the basis function width parameter. An example of a one-dimensional field decomposition is given in Fig. 2. The width and positioning of the basis functions can be determined by spectral analysis (explained in detail in the following section). 8 Amplitude 2 1 0 −1 −10 0 10 Space Figure 2: Example of a one-dimensional field decomposition. A continuous, one-dimensional field decomposed by a finite number of basis functions scaled by the state vector. The dashed line depicts the field, the black solid lines show the Gaussian basis functions and the vertical lines show the position and amplitude of the states. The connectivity kernel can also be decomposed as w (r, r′ ) = ψ ⊤ (r, r′ ) θ, (17) where ψ(r, r′ ) is a vector of Gaussian basis functions and θ is a vector of scaling parameters. A graphical demonstration of the connectivity kernel decomposition and the net shape is provided in Fig. 1. By exploiting the isotropy of the connectivity basis functions, they can be written as ψ(r − r′ ). We will assume that we know the parametric form of the connectivity basis functions, where the scaling parameters θ are unknown. Note, that although most examples presented in this paper use the Mexican hat connectivity kernel, the framework is also valid for other connectivity kernels, provided they are homogeneous and that they can be approximated by a set of isotropic basis functions parametrized by a coefficient vector θ. The mathematical formulation does not restrict the shape of the kernel, where anisotropic arbitrary shapes can be represented. Furthermore, the kernel basis functions need not be a Gaussian shape, but can be, for example, Laplacian to describe the “wizard hat” connectivity or B-spline functions to describe a kernel with a compact support. Each connectivity basis function of the Mexican hat can, individually, be considered a layer in the Wilson and Cowan model. Rearranging the order of the spatial convolution in equation 12 and substituting in equations 15 and 17 to the right hand side we get the approximation of the dynamic field vt+1 (r) ≈ Ts Z f (φ⊤ (r′ )xt )ψ ⊤ (r − r′ ) dr′ θ Ω + ξφ⊤ (r)xt + et (r). 9 (18) Next we cross-multiply equation 18 by φ(r) and integrate over the spatial domain, Ω, to get Z Ω φ (r) vt+1 (r) dr Z Z ψ ⊤ (r − r′ )f (φ⊤ (r′ )xt ) dr′ drθ φ(r) ≈ Ts Ω Ω Z Z +ξ φ(r)φ⊤ (r)drxt + φ(r)et (r) dr. Ω (19) Ω The solution of this equation for xt is equivalent to finding a Galerkin projection, where the field basis functions would be considered the projection functions and the states would be the amplitude coefficients. See Appendix C for more information. In the next step of forming the state-space model we substitute the field decomposition into the left hand side of equation 19 giving Z Ω φ (r) φ⊤ (r) xt+1 dr Z Z ψ ⊤ (r − r′ )f (φ⊤ (r′ )xt ) dr′ drθ φ(r) = Ts Ω Ω Z Z +ξ φ(r)φ⊤ (r)drxt + φ(r)et (r) dr. Ω (20) Ω Now defining the matrix Γ, Z φ (r) φ⊤ (r) dr, (21) Ω and substituting it into equation 20 and cross-multiplying by Γ−1 gives Z ψ ⊤ (r − r′ )f (φ⊤ (r′ )xt )dr′ drθ φ(r) Ω Ω Z −1 φ(r)et (r) dr. + ξxt + Γ xt+1 = Ts Γ−1 Z (22) Ω The analytic derivation of the inner product of n-dimensional Gaussians, which was used for calculating Γ, is provided in Appendix D. Alternative field basis functions may be used to form the decomposition; if the basis is orthonormal, such as a Fourier basis, then Γ is simply the identity matrix. Equation 22 can be simplified by exploiting the isotropy of the connectivity kernel basis functions where ψi (r − r′ ) = ψi (2ci + r′ − r) (23) and ci is the center of the ith basis function of the kernel decomposition. To make the simplification we define [Ψ(r′ )]:i , Ts Γ−1 Z φ(r)ψi (2ci + r′ − r)dr, Ω 10 (24) where [Ψ(r′ )]:i denotes the ith column of Ψ(r′ ) which is a nx × nθ matrix, where nx is the number of basis functions (states) and nθ is the number of connectivity kernel basis functions. This matrix is pre-defined analytically (see Appendix E for analytic convolution of two Gaussians), since it is constant with respect to the states and the unknown parameters. Substituting Ψ(r′ ) into equation 22 gives Z Ψ(r′ )f (φ⊤ (r′ )xt ) dr′ θ Ω Z φ(r)et (r) dr. + ξxt + Γ−1 xt+1 = (25) Ω This substitution swaps the convolution in equation 18, which is computationally demanding, with the convolution in equation 24, which is pre-defined analytically. This is of great importance as it lowers the computational demands of the estimation algorithm and provides a dramatic increase in estimation speed. Now we define the state disturbance as et , Γ−1 Z φ(r)et (r) dr, (26) Ω which is a zero mean normally distributed white noise process with covariance (see Appendix F for the derivation) Σe = Γ−1 Z Z Ω φ (r) γ (r − r′ ) φ (r′ ) ⊤ dr′ drΓ−⊤ . (27) Ω The observation equation of the reduced model is found by substituting equation 15 into equation 13 giving yt = Z m (rn − r′ ) φ⊤ (r′ ) xt dr′ + εt . (28) Ω Note, the sensor index rn has been omitted to delineate the reduced model observation equation from the full model. The observation equation is linear and can be written in the more compact form yt = Cxt + εt , (29) where each element of the observation matrix is Cij , Z Ω m(ri − r′ )φj (r′ ) dr′ . (30) Now we have the final form of the state-space model where xt+1 = Q(xt ) + et (31) yt = Cxt + εt (32) 11 and Q(xt ) = Z Ψ(r′ )f (φ⊤ (r′ )xt ) dr′ θ + ξxt . (33) Ω 2.4. Spectral Analysis Spectral analysis has been used to identify both the number of sensors and the number of basis functions required to reconstruct the neural field from sampled observations (Sanner and Slotine, 1992; Scerri et al., 2009). Based on a two-dimensional extension of Shannon’s sampling theorem (Peterson and Middleton, 1962), the spatial bandwidth of the observed field can be used to provide a lower bound on both the number of sensors and the number of basis functions required to capture the dominant spatial spectral characteristics of the neural field. Let the spectral representation of the post-synaptic membrane voltage field at time t be denoted by Vt (ν), where ν is the spatial frequency. According to Shannon’s sampling theorem, the field must be spatially band-limited for an accurate reconstruction using spatially discrete observations. However, an approximate reconstruction can be obtained if the field is approximately band-limited with Vt (ν) ≈ 0 ∀ν > ν c , (34) where ν c is a cutoff frequency (typically taken as the -3 dB point). Given such a band-limited field, the distance between adjacent sensors, ∆y , must satisfy ∆y ≤ 1 , 2ρy ν c (35) where ρy ∈ R ≥ 1 is an oversampling parameter. A derivation of the sampling theorem that is based on Sanner and Slotine (1992) and Scerri et al. (2009) is provided in Appendix G. The condition in equation 35 must be satisfied to avoid spatial spectral aliasing effects when reconstructing the hidden dynamic field, vt (r), using the sampled observations, yt . In practice, it is difficult to estimate the bandwidth of the cortex using traditional electrophysiological measurements, possibly preventing placement of sensors in accordance to equation 35. However, we envisage it may be possible to estimate the spatial bandwidth using other modalities with higher spatial resolution such as fMRI, NIRS or other optical imaging techniques (Issa et al. (2000)). The expected bandwidth of the human neocortex is approximately 0.1-0.4 cycles/mm (Freeman et al., 2000). Spectral aliasing can be avoided by a proper choice of the spatial sampling distance given the sensors’ spectral characteristics. The spatial extent of the sensors results in a spectral low-pass action, thus providing spatial anti-aliasing filtering. Such sensors limit the spatial bandwidth of the observation to a value denoted by ν cy . Aliasing can then be avoided by positioning the sensors according to equation 35, replacing ν c with ν cy . 12 Although such sensors avoid errors due to aliasing, they attenuate the high spatial frequency variations in the observations. Therefore, any procedure applied to estimate the original field or the underlying connectivity structure from these band-limited observations has the potential to underestimate the high spatial frequency components of both the field and the connectivity kernel. This motivates the use of sensors with wider bandwidths (narrower in space). This choice would require more sensors in order to satisfy Shannon’s sampling theorem for a given spatial region. Therefore, a compromise should be found between the bandwidth of the sensor response, the accuracy of the estimation results, the number of sensors and the computational demands of the estimation procedure when designing experiments. Similar considerations need to be made regarding the representation of the dynamic field, vt (r), using the basis function decomposition. Again, applying Shannon’s sampling theorem, the maximum distance between basis functions to accurately reconstruct the neural field must satisfy 1 2ρφ ν cy ∆φ ≤ (36) where ρφ ∈ R ≥ 1 is an oversampling parameter to determine the basis function separation. The field basis function widths can also be inferred using spectral considerations (Sanner and Slotine, 1992; Scerri et al., 2009). To demonstrate this, we begin with a two-dimensional Gaussian basis function located at the origin, defined as ! 1 ⊤ φ(r) = exp − 2 r r . σφ (37) The corresponding Fourier transform is Φ(ν) =  1 πσν2    1 exp − 2 ν ⊤ ν , σν (38) where σν2 = 1 π 2 σφ2 . (39) To obtain a 3 dB attenuation at ν cy , the basis function width, σφ2 , should be chosen such that σφ2 = where σν2cy = 1 , π 2 σν2cy (40) 2ν ⊤ cy ν cy . ln 2 (41) This ensures that the basis functions can represent the neural field with frequency content up to ν cy . Derivations for equations 38 and 41 are given in Appendix H. Alternatively, the basis function width can be set to provide a desired level of smoothness in the approx13 imated field by reducing the spatial bandwidth. This allows for an approximate representation by a lower order model (fewer basis functions), thus reducing the computational demands of the estimation process. The spatial cutoff frequency of this smoothed field becomes a design choice. To calculate the cutoff frequency imposed by the basis function width, we rearrange equations 40 and 41 giving νcφ 1 = πσφ r ln 2 . 2 (42) Note that for a spatially homogeneous and isotropic field, the basis functions can be placed on a regular grid. Thus, the knowledge of the distance between basis functions can be used to determine the total number of basis functions required to represent a known spatial region. 2.5. State and Parameter Estimation In this section, we describe the procedure for estimating the states, xt , the connectivity kernel parameters, θ, and the synaptic dynamics, ξ. The framework assumes that the parameters of the model (and brain) are stationary over the estimation period. This is a reasonable assumption when applying the framework to an evoked-potential experiment over a short time period in a controlled environment. The estimation process is a two part iterative algorithm, consisting of a state estimation step followed by a parameter estimation step. Once initialized, at each iteration the sequence of estimated state vectors is used to update the parameter set. The resulting parameters are then used when estimating the new state vector sequence for the next iteration. The procedure stops when the parameter estimates converge. The algorithm is initialized using a bounded random state vector sequence, guaranteeing that the initial estimated parameter set forms a stable kernel. The additive form of the unscented Rauch-Tung-Striebel smoother (URTSS) (Sarkka, 2010) was used for the state estimation. The URTSS incorporates an unscented Kalman filter (UKF) (Julier and Uhlmann, 1997; van der Merwe, 2003) in a forward step to estimate posterior states, x̂ft , followed by a backward step to compute the smoothed state estimates, x̂bt . The first and second order moments of the predicted state are captured by propagating the so-called sigma points through the state equation. The sigma points, Xi , are calculated using the unscented transform as follows: X0 = x̄ Xi = x̄ + Xi = x̄ − (43) p p (nx + λ)Px (nx + λ)Px   i , i−nx i = 1, . . . , nx (44) , (45) i = nx + 1, . . . , 2nx where x̄ represents either x̂ft or x̂bt , Px is the corresponding covariance matrix for filtering or smoothing,  p (nx + λ)Px is the ith column of the scaled matrix square root of Px , and nx is the dimension of the i 14 state space. The total number of sigma points is 2nx + 1. The scaling parameter, λ, is defined as λ = α2 (nx + κ) − nx , (46) The positive constant α, that determines the spread of the sigma points around x̄, was set to 10−3 . It is typically set arbitrary small to minimise higher-order effects (Haykin, 2001). The other scaling parameter, κ, was set to 3 − nx in accordance to Julier et al. (2002). The sigma vectors are propagated through the system equations and weighted to form the predicted mean and covariance. The weights are calculated by (m) W0 (c) W0 (m) Wi λ nx + λ λ + (1 − α2 + β) = nx + λ 1 (c) = Wi = i = 1, . . . , 2nx , 2(nx + λ) = (47) (48) (49) where the superscripts m and c denote mean and covariance and β incorporates prior knowledge of the distribution of the states, x (for a Gaussian disturbance, β should be set to 2 (Haykin, 2001)). Since the observation equation is linear (equation 29), the standard Kalman Filter update equations are used to correct the predicted states. The state estimates from forward filtering are used to form a new set of sigma points for the smoother, as described above. The cross-covariance matrix of the states, M, is also required for computing the smoother gain, as described in the summary of the steps in the URTSS algorithm given in Table 2. Although the system is nonlinear, the parameters of the system are linear with respect to the state. This is exploited by our procedure where the parameter estimation uses a least squares (LS) method that minimizes the sum of the squared errors (of a predicted state update) with each new estimate of the state vector sequence. To create the least squares parameter estimator, we first define the nx × nθ matrix q(xt ) = Z Ψ(r′ )f (φ⊤ (r′ )xt ) dr′ . (50) Ω Now given a state estimate sequence from the initialization or an iteration of the URTSS, we can write x̂b1 = q(x̂b0 )θ + ξx̂b0 + e0 x̂b2 = q(x̂b1 )θ + ξx̂b1 + e1 .. . x̂bT = q(x̂bT −1 )θ + ξ x̂bT −1 + eT −1 . 15 1. Forward initialization x̂0 , P0 f 2. Forward iteration: for t ∈ {0, · · · , T }, calculate the sigma points Xi,t using equations 43-46 and propagate through equation 33 f− f Xi,t+1 = Q(Xi,t ) i = 0, . . . , 2nx Calculate the predicted state and the predicted covariance matrix P2n (m) f − − x̂ft+1 = i=0x Wi Xi,t+1 P2n (c) − f− − f− − ⊤ Pft+1 = i=0x Wi (Xi,t+1 − x̂ft+1 )(Xi,t+1 − x̂ft+1 ) + Σe Compute the filter gain, the filtered state and the filtered covariance matrix using the standard Kalman filter update equations − − Kt+1 = Pft+1 C⊤ (CPft+1 C⊤ + Σε )−1 − − x̂ft+1 = x̂ft+1 + Kt+1 (yt+1 − Cx̂ft+1 ) − Pft+1 = (I − Kt+1 C)Pft+1 3. Backward initialization PbT = PfT , x̂bT = x̂fT b 4. Backward iteration: for t ∈ {T − 1, · · · , 0} calculate the sigma points Xi,t and propagate them through equation 33 b− b Xi,t+1 = Q(Xi,t ) i = 0, . . . , 2nx Calculate the predicted state and the predicted covariance matrix P2nx (m) b− x̂b− Xi,t+1 t+1 = i=0 Wi P2nx (c) b− b− b− b− ⊤ Pb− t+1 = i=0 Wi (Xi,t+1 − x̂t+1 )(Xi,t+1 − x̂t+1 ) + Σe P2n (c) b− b ⊤ − x̂ft )(Xi,t+1 − x̂b− Mt+1 = i=0x Wi (Xi,t t+1 ) Compute the smoother gain, the smoothed state and the smoothed covariance matrix −1  St = Mt+1 Pb− t+1   x̂bt = x̂ft + St x̂bt+1 − x̂b− t+1  ⊤  Pbt = Pft + St Pbt+1 − Pb− t+1 St Table 2: Algorithm for the Unscented RTS Smoother. This table shows the steps in the unscented Rauch-Tung-Striebel smoother algorithm. The steps are iterated 10 times for our state estimation procedure. The least squares algorithm is run after each iteration to update the parameter estimates. This can be written in the compact form Z = XW + e, 16 (51) where  x̂b1   x̂b  2 Z=  ..  .  x̂bT      X=       ,    q(x̂b0 ) x̂b0 q(x̂b1 ) .. . x̂b1 .. . q(x̂bT −1 ) x̂bT −1         and   W= θ ξ    e=     , e0 e1 .. . eT −1     .    Following this, the LS parameter estimates (Ljung, 1999), Ŵ, are Ŵ = (X⊤ X)−1 X⊤ Z. (52) 3. Results We now provide demonstrations of the estimation procedure using the new methods that are described above. Data was generated using the neural field model, described in Section 2.2, to demonstrate the state and parameter estimation framework. All parameters for the model are given in Table 3, unless otherwise explicitly stated for particular experiments. Note that the parameters for the spatial and temporal aspects of the model can be considered arbitrary when demonstrating the estimation methodology. This is because the temporal aspects scale with the sampling rate and the length of the times series, and the spatial aspects scale with the size of the neural field and the spatial step size used for numerical simulation. When generating the simulated data, we assumed free boundary conditions for the neural field. An example of the simulated field at a single point in time can be seen in Fig. 3(A). Fig. 3(B) shows an example of data from the modeled sensors and Fig. 3(C) shows the power spectral density (PSD) of the simulated time-series. 2 The sensors were spaced on a 14 × 14 square grid and the sensor width was set to σm = 0.81. This provides a width of 1.5 mm at half the maximum amplitude of the observation kernel. The distance between sensors was set to ∆y = 1.5 mm. This arrangement modeled the recording as having some cross talk (from volume conduction), which is typical for neurophysiological recordings. 3.1. Spatial Frequency Analysis The spectral analysis was used to confirm that the spacing and width of the sensors was adequate to capture the significant dynamics of the field. Fig. 4(A) shows the spatial frequency of the neural field and Fig. 4(B) shows the spatial frequency of the reconstructed neural field from the basis function decomposition. The cutoff frequency (−3 dB point) of the neural field is ν c ≈ 0.24 cycles/mm. From equation 35, this yielded 17 Symbol Quantity Value Units 0.5 0.001 400 mm s n.a. Domain and indices ∆ Ts T τ ς v0 θ σ ψ 1 , σψ 2 , σ ψ 3 ny ∆y σm Σε σγ σd nx ∆φ σφ ρy ρφ α β κ λ spatial discretisation step time step number of time steps used in estimation Model synaptic time constant slope of sigmoidal activation function firing threshold vector of connectivity kernel parameters connectivity kernel width parameters number of sensors distance between adjacent sensors observation kernel width parameter observation noise variance disturbance spatial covariance width disturbance variance Reduced Model number of basis functions distance between field basis functions width of field basis functions in spatial domain Frequency Analysis sensor over-sampling parameter field basis function over-sampling parameter Estimation determines the spread of sigma points prior knowledge of the sigma points distribution scaling parameter (3 − nx ) scaling parameter 0.01 0.56 (Wendling et al., 2005) 1.8 (Marreiros et al., 2008)  ⊤ 100 −80 5 1.8, 2.4, 6 196 1.5 0.9 0.1 × Iny 1.3 0.1 s−1 mV−1 mV mV.spike−1 mm n.a. mm mm mm2 mm mV 81 2.5 1.58 n.a. mm mm2 1 1.67 n.a. n.a. 0.001 2 -78 -80.99 n.a. n.a. n.a. n.a. Table 3: Parameters. Parameters for the neural field model, the reduced model, the frequency analysis and the estimation procedure. 18 A 10 Space 2.4 0 0.0 −1.8 -10 0 B 10 8mv Channels Space C 50ms Power (dB) 45 40 35 30 25 20 15 100 101 102 Frequency (Hz) Figure 3: Example of the spatiotemporal properties of the model and synthetically generated data. (A) The neural field at a single time instant. (B) Data generated by the full neural field model at the first 5 sensors. (C) Mean power spectral density of the data generated by the observations. A 20 30 40 B 30 40 Spatial Freq. 2 Spatial Freq. 2 20 0 Spatial Freq. 0 2 Spatial Freq. 2 Figure 4: Spatial frequency analysis: the neural field. (A) The average (over time) power in dB of the spatial frequency of the neural field. (B) The average (over time) power in dB of the spatial frequency of the reconstructed neural field from the basis function decomposition. The reconstructed field shows a narrower bandwidth as a result of the band limiting properties of the sensors and basis functions. 19 a maximum sensor separation of 2.08 mm. This confirms that the separation of ∆y = 1.5 mm was sufficient to mitigate problems associated with spatial aliasing of the dominant dynamics. The spatial frequency of the observed field is shown in Fig. 5(A). The cutoff frequency is ν cy ≈ 0.2 cy- A 29 B C 20 20 νcy ≈ 0.2 0.0 0.6 Spatial Freq. 23 νcφ ≈ 0.12 Power (dB) 27 Power (dB) Spatial Freq. 0.6 −20 −60 −100 0.0 0.5 Spatial Freq. 1.0 −20 −60 −100 0.0 0.5 Spatial Freq. Figure 5: Spatial frequency analysis: observations and basis functions. The dashed line shows the cutoff frequency (−3 dB point). (A) The average (over time) power in dB of the spatial frequency of the observations. (B) The spatial frequency response of a one-dimensional model sensor. The symmetry of the sensors and field basis functions allow for the one-dimensional representation of the spectral characteristics. (C) The spatial frequency response of a one-dimensional neural field basis function, illustrating a narrower band width than the sensors. cles/mm (−3 dB point). The observation cutoff frequency is lower than the field cutoff due to the low-pass filtering effect of the sensors. The spatial frequency response of the sensors can be seen in Fig. 5(B). As described in Section 2.4, the basis function configuration can be systematically chosen to account for the full spatial bandwidth from the observations, or alternatively, the configuration can be chosen to represent the neural field up to a specified bandwidth to limit the computational complexity of the estimation procedure. In the examples provided in this paper, the later option was chosen, where the neural field cutoff frequency was further reduced to ν cφ = 0.12 cycles/mm. The band-limiting effect on the field reconstruction can be seen by comparing Fig. 4(A) to Fig. 4(B). This effect is due to the spatial frequency response of the basis functions, which is shown in Fig. 5(C). Due to the slow roll-off in the frequency response of Gaussian basis functions an oversampling parameter of ρφ = 1.67 was chosen. Given the spectral response of the basis functions, the over-sampling parameter, and the spatial domain an equally spaced grid of 9 × 9 basis functions was used to decompose the field, which resulted in 81 states. 3.2. State and Parameter Estimates Several experiments were performed to demonstrate the performance of the state and parameter estimation algorithm. The first experiment evaluated the estimation performance using a Monte Carlo approach for a single set of parameters. In this experiment we demonstrate convergence of the algorithm, show the distribution of estimated parameters and the accuracy of the state estimates for a single set of parameters. 20 1.0 The second experiment evaluated the performance with a selection of different connectivity kernel parameters. In addition, results are shown where the support of the connectivity kernels differs between actual data and the estimator. The third experiment shows how the estimation accuracy depends on the parameters of the sigmoidal nonlinearity. The final experiment shows the effect of the width of the observation kernel on the estimation performance. 3.2.1. Experiment I: Monte Carlo Simulation with Fixed Parameters A Monte Carlo approach was employed where 150 realizations of the data were generated, such that distributions of the parameters estimates could be formed. Each realization consisted of 500 ms of data (sampled at 1 kHz) and the estimation was applied to the final 400 ms, allowing the model’s initial transients to die out. The initial state and parameters were unknown to the estimator. Ten iterations of the estimation algorithm were used to ensure that the parameter estimates converged. Fig. 6 shows the rate of convergence, confirming that the algorithm converged to steady parameter values, with the changes in the parameter errors (between consecutive iterations) falling to less than 10−4 after 6 iterations. 30 ξ θ0 θ1 θ2 25 Error 20 15 10 5 0 1 2 3 4 5 6 7 8 9 10 Iteration Figure 6: Convergence of parameters. The mean of the absolute error between true and estimated parameters averaged over 150 realizations. All parameters converge after 6 iterations, where the mean difference between errors in consecutive iterations falls below 10−4 . Histograms of the parameter estimates obtained over the 150 realizations are shown in Fig. 7. In each of the subplots of the figure the true parameter value is shown by the solid line and the mean parameter estimate is shown by the dashed line. The actual parameter values, the mean parameter estimates, standard deviations, and biases for this parameter combination are given as the first entry in Table 4, where other entries in the table show results from Experiment II. The biases show the percentage differences between the actual values and the mean parameter estimates. For this parameter combination the mean estimates are in good agreement (within 1 standard deviation) with the actual values of the connectivity kernel parameters. However, the time constant parameter, ξ, is outside the estimated distribution (with bias of 2.67%), with a small variance in parameter estimates. We attribute this small bias to the loss of the high frequency 21 A B 40 40 20 20 0 40 100 C θ̂0 160 0 −130 −90 D θ̂1 0 0.85 0.90 40 70 20 35 0 3.0 5.0 6.5 −50 0.95 ξˆ θ̂2 Figure 7: Histograms of the parameter estimates over 150 realizations. In each of the subplots the solid lines show the actual parameter values and the dotted lines show the means of the estimated parameters. (A) The histogram of the central excitatory connectivity kernel basis function parameter estimates, θ0 . (B) The histogram of the surround inhibition connectivity kernel basis function parameter estimates, θ1 . (C) The histogram of the longer range excitatory connectivity kernel basis function parameter estimates, θ2 . (D) The histogram of the ξ estimates where ξ = 1 − ζTs and ζ is the reciprocal of the synaptic time constant. Experiment I and II-A II-B II-C Parameter θ0 θ1 θ2 ξ θ0 θ1 θ2 ξ θ0 θ1 θ2 ξ Actual 100 -80 5 0.9 90 -80 5 0.9 110 -80 5 0.9 Mean Est. 101.75 -81.00 4.76 0.924 103.86 -89.05 5.23 0.923 98.90 −73.13 4.37 0.925 Std. Dev. 21.30 14.82 0.65 0.003 23.14 16.57 0.74 0.003 16.48 11.72 0.54 0.003 Bias (%) 1.75 1.25 4.8 2.67 15.4 11.31 4.6 2.56 10.1 8.58 12.6 2.78 Table 4: Results from varying the Mexican hat connectivity kernel parameters 22 components in the smoother estimated field with respect to the actual field. These results indicate that choosing to limit the frequency content of the reconstructed neural field has only minor effects on the accuracy of the estimated kernel while allowing for significant gains in the computational demand. The accuracy of the state estimates was evaluated by comparing the field reconstruction to the true field using the mean (over iterations) of the root mean square error (MRMSE) over space. The MRMSE of the field estimates over time for 150 realizations is shown in Fig. 8 with the 95 % confidence interval. At each time point, the error is consistently around 0.5 mV, which is approximately 10.73 % of the field 0.52 MRMSE 0.51 0.50 0.49 0.1 0.2 0.3 0.4 0.5 Time (sec) Figure 8: Error in the field reconstruction. The mean (over iterations) of average (over space) RMSE of the estimated field over 150 realizations (solid line) and the 95 % confidence interval (grey region). voltage range ([-2.20, 2.46] mV). The error in the reconstructed field is again mostly due to the loss of high spatial frequency information from the basis function decomposition as illustrated in Fig. 9. The figure C B 10 2.0 0 0.0 Space Space 10 0 0.0 −1.5 −10 0 Space 10 10 2.0 Space A 2.0 0 0.0 −1.5 −10 0 Space 10 −1.5 −10 0 10 Space Figure 9: Examples of the actual and estimated neural fields and the error. (A) The actual field at a single time instant from the neural field model that was used to generate the data. (B) The reconstructed field of the reduced model, showing the effect of the basis function decomposition on the frequency content of the estimated field. (C) The error between the true and reconstructed fields. shows a snapshot at a single time instant of the true neural field, the estimated field, and the corresponding error, which is the difference between the true and estimated neural fields. It can be seen in Fig. 9(B) that 23 the estimated field contains the basic structure of the true field, shown in Fig. 9(A), but is missing the higher-frequency detail, which can be see in the error shown in Fig. 9(C). 3.2.2. Experiment II: Results for Different Connectivity Kernel Parameters To study the robustness of the estimation algorithm, the performance was evaluated for various connectivity kernel parameters. The reconstructed connectivity kernels from the estimated parameters are shown in Fig. 10. For Fig. 10(A-C) Mexican hat connectivity kernels were used. For each subset of parameters, Connection Strength A 40 θ0 =100 θ1 =-80 θ2 =5 θ0 =90 θ1 =-80 θ2 =5 θ0 =110 θ1 =-80 θ2 =5 0 −20 D Connection Strength C B E 60 F θ0 =15 θ1 =0 θ2 =0 θ0 =30 θ1 =0 θ2 =0 θ0 =50 θ1 =0 θ2 =0 0 0 −10 10 0 −10 Space Space 10 −10 0 Space Figure 10: Cross-section of connectivity kernel estimates. In each of the subplots the actual kernels are shown by the solid lines, the mean kernel estimates (over 50 realizations) are shown by the dotted lines and the 95 % confidence intervals are shown by the shaded grey regions. the mean of the reconstructed kernel over 50 realizations is shown as a dotted line and the 95 % confidence interval is indicated by the grey shaded region. In each case, the actual kernel was within the 95 % confidence interval of the kernel estimate. The statistics for the kernel parameter estimates for each of the parameter sets are shown in Table 4. For each of the Mexican hat connectivities, the support of the kernel was assumed to be known. However, in Fig. 10(D-F) the kernel support in the estimator differed from the actual data. In these cases the data sets were generated using a single positive (excitatory) Gaussian kernel, as seen by the solid line in each of the subplots of the figure. The estimation algorithm was constructed allowing for the Mexican hat connectivity of the previous examples. Each plot shows that the actual kernel lies within the confidence intervals of the reconstructed estimates for the majority of the kernel space, until the amplitude is relatively low. 24 10 Table 5 shows the statistics of the parameter estimates. Although the net reconstructed kernel estimates Experiment II-D II-E II-F Parameter θ0 θ1 θ2 ξ θ0 θ1 θ2 ξ θ0 θ1 θ2 ξ Actual 15 0 0 0.9 30 0 0 0.9 50 0 0 0.9 Mean Est. 35.68 −18.16 0.81 0.923 40.43 −12.01 0.35 0.924 61.46 −15.03 0.27 0.924 Std. Dev. 16.72 11.77 0.53 0.002 11.47 7.8 0.36 0.003 11.51 7.56 0.29 0.003 Bias (%) 137.87 2.56 34.77 2.67 22.92 2.67 Table 5: Results from varying the Gaussian connectivity kernel parameter are in good accordance with the actual kernels, the individual kernel basis function weights are not in agreement with the actual parameters. This result demonstrates that the interpretation of the results should be based on the net shape of the estimated connectivity kernel and not on the estimated parameter values individually. These results indicate that the solution to the system may not be unique in non-ideal, practical situations. Nevertheless, the net shape of the connectivity kernel can be determined reasonably accurately. 3.2.3. Experiment III: Results for Varying the Parameters of the Sigmoidal Activation Function Fig. 11 demonstrates the effect of the sigmoid parameters on the estimation accuracy. The RMSE of the reconstructed connectivity kernel was calculated for various parameters of the sigmoid by varying the slope, ς, and the threshold, v0 . The other model parameters remained fixed and were set according to the Table 3. From Fig. 11 it can be seen that the estimation accuracy depends on the parameters of the sigmoid function. This is expected since that by altering the parameters we are changing the nonlinearity and the saturation region of the sigmoid. For example, when the sigmoid slope is very gentle, most of the neural activity falls in the central, nearly linear region of the sigmoid, thus having minor effects on the estimation procedure. On the other hand, sharp slopes result in more of the activity falling in the saturation region of the sigmoid with adverse effects on the estimation. Similar arguments can be made for variations to the sigmoid threshold. 3.2.4. Experiment IV: Results for Varying the Sensor Kernel Width The accuracy of the estimation algorithm was tested for different observation kernels by varying the width parameter, σm . Fig. 12 shows the RMSE of the connectivity kernel estimates for the various observation kernel widths. As the theory presented earlier in Section 2.4 states, the observation kernel acts as a spatial 25 Figure 11: RMSE of connectivity kernel estimates with various sigmoid parameters. Log of the RMSE of the connectivity kernel for different values of firing threshold, v0 , and the slope of the sigmoid, ς. Mean kernel estimates are obtained over 50 realizations in each case. White spaces indicate cases where the kernel can not be estimated due to numerical problems. 6 5 RMSE 4 3 2 1 0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0 3.3 3.6 3.9 4.2 4.5 Sensor Kernel Width Figure 12: RMSE of connectivity kernel estimates with various observation kernel widths. RMSE of the connectivity kernel for different values of observation width, σm . Mean kernel estimates are obtained over 50 realizations in each case. 26 filter. When the observation kernel is too narrow, high frequency components in the neural field are present but under sampled. These high frequency components result in aliasing errors in the estimated field, thus affecting the kernel estimation accuracy. On the other hand, for the very wide observation kernels considered, important neural dynamics are not represented in the estimated field due to the filtering out of the important frequency components by the sensors, resulting in poor accuracy in the kernel estimate. As indicated in Fig. 12, the best results are obtained when the sensor width is chosen to minimize both sources of errors. 4. Discussion 4.1. Implications for Experimental Design As seen in Fig. 7 and Fig. 10, the results show good estimation accuracy for the connectivity kernels, where the true kernels lie within the 95 % confidence intervals of the kernel estimates. The estimates of the connectivity kernels show relatively broad distributions over the different realizations, which can be expected considering the system’s input is modeled as a random disturbance and the observations are corrupted by noise. To overcome the effect of noise, the use of a controlled input stimulus should be considered when applying the estimation algorithm to electrophysiological recordings, where given knowledge of the input and the appropriate model, it is expected that the variance of the parameter estimates would decrease. Even when using a known input, our results imply that we should only expect to learn distributions of parameters, rather than uncover a ‘true’ parameter set. By designing an input signal carefully and using an evoked-potential paradigm where stimuli are repeated many times, recordings could be averaged over a sufficient number of trials (realizations), such that a distribution of the parameter estimates and resulting connectivity kernels can be formed. The frequency analysis specifies the minimum sensor arrangement that is required to capture the dominant cortical dynamics. As pointed out in numerous publications by Freeman and colleagues (see Freeman and Baird (1987) for example), a detailed description of the spatial properties of responses to stimuli is essential in understanding cortical dynamics. Importantly, the spatial frequency analysis tools presented in this paper can be used to design high density intracranial electrophysiological recording systems to avoid spatial aliasing, ensuring the essential spatial dynamics can be observed. In addition, the frequency analysis highlights the relationship between the width of the sensors and the information in the data. This analysis is becoming increasingly relevant as recording systems are becoming more sophisticated with higher density electrodes (Brinkmann et al., 2009). The work of Freeman et al. (2000) provides direct evidence that the spatial frequency of the neural field is in the range of 0.1-0.4 cycles/mm and spacing required to prevent aliasing is approximately 1.25 mm. Micro-electrodes systems that satisfy this requirement are now currently in use for brain computer interface 27 and epilepsy monitoring applications (Cheung, 2007), which suggests that the framework can be applied to human data with existing technology. 4.2. Computational Complexity The closest competitor to our ‘grey-box’ integro-difference equation approach is perhaps the multivariate auto-regressive (MVAR) model, which is a ‘black-box’ approach. Examples of MVAR models for estimating functional connectivity include Granger causality (Hesse et al., 2003), the direct transfer function (Kaminski and Blinowska, 1991), and partial directed coherence (Sameshima and Baccalá, 1999). A distinct advantage of the integro-difference equation based technique proposed in this current study is that the number of parameters required to model the system is considerably less than a typical MVAR model. Furthermore, by using the neural field equations as a parametric form of a model, the estimated quantities have physical meaning and, therefore, can be clinically relevant. An advantage of the MVAR-based methods is that they allow for inhomogeneous connectivity, which is present in the longer range cortico-cortical fibres. Previous work by Schiff and Sauer (2008) demonstrated for the first time the applicability of the unscented Kalman filter (UKF) for state estimation of discrete (in space and time) nonlinear spatiotemporal systems. Within their framework, each discrete spatial location formed an element of the state vector to be estimated with the UKF. A limitation of this framework was that an increase in spatial resolution and/or the physical size of the neural field resulted in an increased state dimension. This, in turn, leads to increases in uncertainty in the state estimates. In addition, directly linking the complexity of the underlying system to the spatial resolution of the observation process implies that the system’s dynamics become increasingly more complex the closer we look, which may not be an accurate modeling assumption. The current study builds on Schiff and Sauer’s approach by estimating a continuous field whose complexity is independent of the resolution of the observation process and spatial resolution of the model. Instead, the dimensionality is governed by the expected complexity of the underlying field, which we measure using the frequency analysis presented in the Methods section. These methods provide a rigorous framework for determining the dimensionality of the system using inherent spatial properties of the field. It should be noted that Schiff and Sauer suggest the possibility of using a Galerkin projection, or similar decomposition that is used in this paper, to perform efficient estimation as a direction for future work. Indeed, we have shown that this enables estimation of a larger field with a lower dimensional state vector than would be otherwise possible. The isotropy of the connectivity kernel basis functions means that the mixing effect of the kernel on the firing rates can be written as a convolution of the kernel with the field basis functions. This transformation greatly simplifies the computational complexity of the estimator, which is important for the estimator to run in a real-time scenario. 28 4.3. Patient-Specific Modeling Generating data-driven neural field models has the potential to have significant impact on several areas of neuroscientific research and clinical neurology. Specifically, these subject-specific neural field models will provide new avenues for the application of engineering techniques to neural systems. For example, state tracking could be used in brain-computer interfaces or epileptic seizure prediction, where specific regions of the state-space or certain parameter combinations may indicate intent of movement or imminent seizures, respectively. In addition, a patient-specific state-space model will allow application of electrical stimulation to robustly prevent the neurodynamics entering pathological regions of state-space, thereby possibly preventing seizure onset or abating seizure activity. Estimation of the connectivity structure or the synaptic time constant may also have significant impact on the treatment of disease. For example, recent theoretical work has demonstrated that seemingly similar electrographic seizures in absence epilepsy may arise from either abnormal excitatory or inhibitory mechanisms (Marten et al., 2009). These very different mechanisms would require very different medications for successful treatment. However, the underlying mechanisms would be hidden from the clinician in the normal clinical setting. Successful parameter estimation has the potential to reveal these hidden mechanisms allowing for improved treatment strategies. The development of the estimation framework also has implications for theoretical neuroscience, where the ability to estimate subject-specific properties will allow testing of hypotheses that have been generated in computational modeling studies. This has the potential to validate neural field models and, in doing so, provide new insights that may lead to a greater understanding of mesoscopic and macroscopic neurodynamics. For example, using the proposed model-based framework, we can reconstruct a neural field from data and estimate parameters to check if they correspond to theoretically derived values where certain solutions may exist. In addition, the parameters and the state-space can be explored during cognitive tasks and states of arousal such as sleep and anaesthesia. 4.4. Extensions to the Framework In this study, an assumption was made that the chosen kernel basis functions can exactly represent the actual connectivity kernel. In future work, we plan to relax this assumption and show how the support of the connectivity kernel can also be inferred from data. In addition, we have assumed that the neural field is isotropic and homogeneous. It is generally accepted that the intracortical connectivity structure can be accurately modeled as being isotropic. However, the assumption of homogeneity is quite strong as the size of the neural field becomes larger. Future work should focus on developing estimation algorithms where assumptions of homogeneity can be relaxed. Other natural extensions of this framework are the inclusion of a post-synaptic response kernel with a finite rise and decay time, and a finite distant-dependant action potential propagation velocity into the parametric neural field equations. 29 The assumption was made that the parameters of the model are stationary over the estimation period. However, it is well known that factors that influence excitability are non-stationary over longer time periods (see Ullah and Schiff (2010)). Therefore, the estimation framework must be altered if it is to be applied to data over extended periods of time. A possible way to achieve this is to use a more sophisticated model that accounts for longer-term time dependancies (higher temporal order). Alternatively, the parameters of the model could be considered time-varying, where changes in excitability would result in different connectivity kernel parameters. Furthermore, the model may be considered heteroscedastic, where the disturbance is also time-varying. This is common in spatio-temporal econometric studies, where methods have been developed to test for heteroscedasticity (Anselin, 1988). Following this, the disturbance properties of the system will also need to be estimated. A possible approach allowing for time-varying parameters is to perform a windowbased or fixed-lag estimation, which would require the assumption that the parameters are quasi-stationary over the duration of the window. The estimation framework that we have presented is inherently a batch algorithm designed to operate after data has been collected. A further extension of the framework is to allow for the algorithm to learn online. A filter that is designed to be implemented in an online setting can be used in place of the smoother along with an online parameter estimator to update the parameters after a new observation arrived. A real-time implementation can be achieved using the UKF (instead of the URTSS) and an online recursive LS estimator Ljung (1999). 5. Conclusion This paper has presented a framework for creating data-driven neural field models by developing a method for estimating intracortical connectivity and synaptic dynamics from electrophysiological data. In doing so, we have highlighted important aspects of the experiment design process, demonstrating the relationship between sensor bandwidth and the complexity of the field that can be inferred. We have derived a finiteorder state-space representation of the continuum neural field equations that allows inference to be performed efficiently. Importantly, the mathematical formulation that we have presented scales well with the number of sensors. The major contribution of this framework is to formally bridge a gap between experimental and theoretical neuroscience in order to gain a better understanding of the mesoscopic dynamics of the brain. A validation of the estimation framework on real data is now required. 6. Acknowledgments Dean Freestone would like to thank Dr. Mark van Rossum for hosting him while visiting the University of Edinburgh where this work was undertaken and Dr. Matthias Hennig for useful discussions and feedback. He would also like to thank Dr. James Fallon for assistance wit the observation kernel experiment. This 30 research was funded in part by the Australian Research Council (Linkage Project LP0560684). The Bionic Ear Institute acknowledges the support it receives from the Victorian State Government through the Operational Infrastructure Support Program. Dean Freestone would like to thank the University of Melbourne’s Scholarships Office for the Postgraduate Overseas Research Experience Scholarship and the Harold Mitchell Foundation for a traveling scholarship for supporting this research. Parham Aram would also like to thank Andrew Hills for useful discussions and assistance with code and the generation of figures. The authors would also like to thank Dr. John Terry for insightful discussions and feedback. Finally, the authors would like to thank the anonymous reviewers who have provided a critique that has led to valuable insights that would have gone otherwise astray. Appendix A. Time Discretization of the Neural Field Model This appendix provides a description of the discretization of the neural field model to form the integrodifference equation (IDE). To perform the discretization, we used a one-step Euler method, where equation 11 can be approximated by v (r, t + Ts ) − v (r, t) = −ζv (r, t) Ts Z w (r − r′ ) f (v (r′ , t)) dr′ . + (A.1) Ω For clarity, we shall index time points in the discrete time form of the model using the subscript t and the next time point as t + 1. Rearranging equation A.1, we get vt+1 (r) = vt (r) − Ts ζvt (r) + Ts Z w (r − r′ ) f (vt (r′ )) dr′ . (A.2) Ω Now making simplifications, the discrete time form of the model is vt+1 (r) = ξvt (r) + Ts Z w (r − r′ ) f (vt (r′ )) dr′ , (A.3) Ω where ξ = 1 − Ts ζ. Note, the discretization using the first order Euler method has the potential to affect the system’s behaviour, where the difference between the continuous and discrete-time systems can potentially be significant. The sampling period of 1 ms was used in this paper which is ten times smaller than the synaptic time constant, therefore the discrepancy between the continuous and its discrete approximation is small. This yielded a stable discrete-time system and good estimation results. An analysis of the error introduced by Euler discretization method can be found in Butcher and Corporation (2008). An alternative to Euler’s method is the local linearization method (Ozaki, 1994). 31 Appendix B. Observation Kernel The Gaussian shape of the observation kernel was verified experimentally. This was achieved by moving a dipole across a recording electrode submerged in saline. A photograph through a microscope of the submerged electrodes can be seen in Fig. B.13. The measurement sensor was the most central electrode of Figure B.13: Observation kernel experimental setup. Photograph through a microscope of the Utah array recording electrode and the twisted bipolar current source. a 5 × 5 Utah array (Blackrock Microsystems) which had electrode spacing of 400 µm. The current source was a twisted wire electrode (Plastics One Inc.) with tip-tip spacing of 60 µm and tip diameter of 200 µm. The twisted wire electrode was moved in the saline solution using a 3 dimensional manipulator (Sutter Instrument Company, MPC-100) in increments of 200 µm. The twisted wire electrode and the Utah array were co-registered visually using the microscope. The stimulation waveform for the current source consisted of a 40 Hz sine wave with a duration of 1 s. The recording from the central electrode in the array was used to estimate the shape of the sensor kernel. To extract the amplitude of the signal, the data was first band-pass filtered with a center frequency of 40 Hz. The envelope of the signal was found by computing the mean instantaneous amplitude over the 1 s stimulation duration by a= q 1X (yt2 − Hyt2 ) T t (B.1) where H denotes the Hilbert transform operator, t indexes samples and T is the total number of samples in the time series. Fig. B.14 shows the mean amplitude of the recorded envelope versus the distance between the twisted wire electrode and the recording electrode. A Gaussian fitted to the amplitudes is also shown. The figure shows that a Gaussian is indeed a good model for the sensor kernel in the homogeneous solution. 32 ï5 3.5 x 10 measured data fitted Gaussian amplitude (V) 3.4 3.3 3.2 3.1 3 ï2000 0 2000 position (µm) Figure B.14: Observation kernel experimental results. Amplitude of the envelope of the signal recorded from the electrode at the center of the array vs. the source position. Appendix C. Equivalence to Galerkin’s Projection To show how equation 19 is equivalent to the Galerkin projection, we define the residual between the field and the decomposition as Ξ(xt , vt+1 (r)) = vt+1 (r) Z f (φ⊤ (r′ )xt )ψ ⊤ (r − r′ ) dr′ θ − Ts Ω ⊤ − ξφ (r)xt − et (r). (C.1) The Galerkin projection finds xt , such that the inner product of the residual with all of the basis functions satisfies Gi = Z Ξ(xt , vt+1 (r))φi (r)dr = 0, (C.2) Ω where i = 0, 1, . . . , nx − 1, and nx is the number of states and basis functions. If equation C.2 is satisfied the basis functions are orthogonal to the residual. By substituting equation C.1 into equation C.2 and rearranging we get equation 19 of the main text. Note, if Galerkin’s method was used to estimate the states, nx coupled ODEs would need to be solved independently (each Gi ) to estimate the state vector, which is a standard method in fluid dynamics and economics modeling. The method in this paper takes a different approach by forming the state-space model, allowing for more efficient estimation. 33 Appendix D. Inner Product of Two n-Dimensional Gaussians The derivation in this appendix was used for calculating the matrix Γ from equation 21. Consider two Gaussian basis functions and   1 ϕi (r) = exp − 2 (r − µi )⊤ (r − µi ) σi (D.1) ! 1 ϕj (r) = exp − 2 (r − µj )⊤ (r − µj ) . σj (D.2) The inner product of the two n-dimensional Gaussians is defined as ϕi ⊗ ϕ j , Z ϕi (r)ϕj (r)dr. (D.3) Rn A Gaussian basis function is symmetric with respect to its center, therefore ϕi (µi + r) = ϕi (µi − r). (D.4) Letting r = r − µi , equation D.4 becomes ϕi (r) = ϕi (2µi − r). (D.5) Now substituting equation D.5 into equation D.3 gives ϕ i ⊗ ϕj = = Z ϕi (r)ϕj (r)dr ZR n ϕi (2µi − r)ϕj (r)dr Rn = (ϕi ∗ ϕj )(2µi ). (D.6) Now we solve the convolution using the formula derived in Appendix E, and evaluating it at r = 2µi gives ϕi ⊗ ϕ j = πσi2 σj2 σi2 + σj2 ! n2 ! 1 ⊤ × exp − 2 (µ − µj ) (µi − µj ) . σi + σj2 i This completes the derivation. 34 (D.7) Appendix E. Convolution of Two n-Dimensional Gaussians In this section, we provide a derivation for the convolution of two n-dimensional Gaussian basis functions. This derivation is used in the analytic calculation of Ψ(r) and the elements of the observation matrix, C. Consider the two Gaussian basis functions from equations D.1 and D.2. The derivation proceeds by assuming the Gaussian basis functions are positioned at the origin (µi = µj = 0). The convolution theorem states that the Fourier transform of the convolution of two functions is equal to the product of their Fourier transforms. Using the formula from the derivation of the Fourier transform of an n-dimensional Gaussian from equation H.4 in Appendix H, we have F (ϕi ∗ ϕj ) (r) = π 2 σi2 σj2
n2
exp −(σi2 + σj2 )π 2 ν ⊤ ν , where F denotes the Fourier transform. Multiplying and dividing by π σi2 + σj2 written as F (ϕi ∗ ϕj ) (r) = πσi2 σj2 1 σi2 + σj2
n2 ! n2 (E.1)

n2 , equation E.1 can be (π(σi2 + σj2 )) 2
× exp −(σi2 + σj2 )π 2 ν ⊤ ν . n (E.2) Taking the inverse Fourier transform, the convolution becomes (ϕi ∗ ϕj ) (r) = πσi2 σj2 σi2 + σj2 ! n2 ! 1 ⊤ r r . exp − 2 σi + σj2 (E.3) In the case where µi 6= µj 6= 0, equation E.3 becomes (ϕi ∗ ϕj ) (r) = πσi2 σj2 σi2 + σj2 ! n2 ! 1 ⊤ ×exp − 2 (r − µ) (r − µ) , σi + σj2 (E.4) where µ = µi + µj . Appendix F. Mean and Covariance of the State Disturbance In this section, we show that the state disturbance vector is et ∼ N (0, Σe ) if et (r) ∼ GP (0, γ(r − r′ )). We begin the derivation by acknowledging the state disturbance vector is a linear function of et (r). Therefore, 35 et also follows a Gaussian distribution. From equation 26, the state disturbance vector is defined as et , Γ−1 Z φ(r)et (r)dr. (F.1) Ω The expected value of et is given by E [et ] = Γ−1 Z φ (r) E [et (r)] dr = 0. (F.2) Ω The covariance of et is   Σ e = E et e⊤ t Z Z −1 = Γ E[ φ(r)et (r)dr φ⊤ (r′ )et (r′ )dr′ ]Γ−⊤ Ω Z ZΩ −1 φ(r)E[et (r)et (r′ )]φ⊤ (r′ )dr′ drΓ−⊤ =Γ ZΩ ZΩ −1 =Γ φ(r)γ(r − r′ )φ⊤ (r′ )dr′ drΓ−⊤ . Ω (F.3) Ω Appendix G. Derivation of Sampling Theorem This appendix provides a derivation for equations 35 and 36, which provide conditions on distance between the sensors and the field basis functions. The derivation is adapted from Sanner and Slotine (1992) and Scerri et al. (2009). Consider the cortical sheet defined on the surface Ω with sensors placed on a regular grid at positions rn (where n = 0, ..., ny − 1 indexes the sensors) that are separated by the distance ∆y . Noise-free observations can be considered samples from a continuous field that is filtered (convolved) with the observation kernel described by y(rn ) = δ(r − rn )vm (r), (G.1) where vm (r) = vt (r) ∗ m(r). The observations are described in matrix form by ny −1 y= X δ(r − rn )vm (r). (G.2) n=0 The Fourier transform of y is  ny −1  1 X kn Y (ν) = 2 ∗ Vm (ν) δ ν− ∆y n=0 ∆y   ny −1 kn 1 X 1 , Vm ν − = 2 Vm (ν) + 2 ∆y ∆y n=1 ∆y 36 (G.3) (G.4) where k is an ny × 2 matrix, spatially indexing the sensors (the columns are different spatial dimensions such that rn = ∆y kn ). It should be clear that Y (ν) contains Vm (ν) plus ny − 1 copies (aliases) of Vm (ν) repeated on the regular lattice, where each copy is centered at kn ∆−1 y . Due to the lowpass action of the observation kernel, the filtered dynamic field is band-limited with a cutoff frequency ν cy . In order for Vm (ν) to be disjointed from its aliases, the amount of shift ∆−1 must be at least 2ν cy , since each alias has a y spectral range of   kn kn ν= − ν cy , + ν cy . ∆y ∆y (G.5) Therefore, to avoid aliasing the distance between observations must satisfy 1 ≥ 2ν cy . ∆y (G.6) Since the observation kernel is not an ideal filter, an oversampling parameter should be included to further separate the aliases. Rearranging equation G.6 and adding the oversampling parameter, ρ, gives ∆y ≤ 1 , 2ρν cy (G.7) where ρ ∈ R ≥ 1. Appendix H. Spatial Cut-off Frequency of the Field Basis Functions This section provides a derivation for the width of the Fourier transform of the field basis functions in terms of their spatial cut-off frequency (equations 38 and 41). We define an n-dimensional Gaussian centered at the origin as ! 1 ⊤ φ(r) = exp − 2 r r . σφ (H.1) Now applying the n-dimensional Fourier transform yields !
1 ⊤ exp − 2 r r exp −2πiν ⊤ r dr Φ(ν) = σ n R φ ! Z 1 ⊤ exp − 2 [r + σφ πiν] [r + σφ πiν] dr = σφ Rn
× exp −σφ2 π 2 ν ⊤ ν . Z (H.2) Evaluating the integral on the RHS of equation H.2 gives Z !
n 1 ⊤ exp − 2 [r + σφ πiν] [r + σφ πiν] dr = πσφ2 2 . σφ Rn 37 (H.3) Substituting equation H.3 into equation H.2 gives the scaled Gaussian Φ(ν) =  1 πσν2  n2   1 exp − 2 ν ⊤ ν , σν (H.4) where 1 σν2 = π 2 σφ2 . (H.5) For 3 dB attenuation (1/2 to peak power) at the cut-off frequency, ν c , we need to set 2 |Φ(ν c )| = 1 2 |Φ(0)| , 2 (H.6) since the maximum of the Gaussian described in equation H.4 is at the origin. Solving for σν2 , we get σν2 = 2ν ⊤ c νc . ln 2 38 (H.7) References Amari, S., 1977. Dynamics of pattern formation in lateral-inhibition type neural fields. Biological Cybernetics 27, 77–87. Anselin, L., 1988. Spatial econometrics . Braitenberg, V., Schüz, A., 1998. Cortex: statistics and geometry of neuronal connectivity. Springer Berlin. Breakspear, M., Roberts, J., Terry, J., Rodrigues, S., Mahant, N., Robinson, P., 2006. A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cerebral Cortex 16, 1296. Bressloff, P., Cowan, J., Golubitsky, M., Thomas, P., Wiener, M., 2001. Geometric visual hallucinations, Euclidean symmetry and the functional architecture of striate cortex. Philosophical Transactions B 356, 299. Brinkmann, B., Bower, M., Stengel, K., Worrell, G., Stead, M., 2009. Large-scale electrophysiology: acquisition, compression, encryption, and storage of big data. Journal of Neuroscience Methods 180, 185–192. Butcher, J., Corporation, E., 2008. Numerical methods for ordinary differential equations. Wiley Online Library. Cheung, K., 2007. Implantable microscale neural interfaces. Biomedical microdevices 9, 923–938. Coombes, S., 2005. Waves, bumps, and patterns in neural field theories. Biological Cybernetics 93, 91–108. Coombes, S., Venkov, N.A., Shiau, L., Bojak, I., Liley, D.T.J., Laing, C.R., 2007. Modeling electrocortical activity through improved local approximations of integral neural field equations. Physical Review E 76, 051901. Daunizeau, J., Kiebel, S., Friston, K., 2009. Dynamic causal modelling of distributed electromagnetic responses. NeuroImage 47, 590–601. David, O., Friston, K., 2003. A neural mass model for MEG/EEG: coupling and neuronal dynamics. NeuroImage 20, 1743–1755. David, O., Kiebel, S., Harrison, L., Mattout, J., Kilner, J., Friston, K., 2006. Dynamic causal modeling of evoked responses in EEG and MEG. NeuroImage 30, 1255–1272. Deco, G., Jirsa, V., Robinson, P., Breakspear, M., Friston, K., 2008. The dynamic brain: From spiking neurons to neural masses and cortical fields. PLoS Computational Biology 4. Dewar, M., Scerri, K., Kadirkamanathan, V., 2009. Data-driven spatio-temporal modeling using the integro-difference equation. IEEE Transactions on Signal Processing 57, 83–91. Freeman, W., Baird, B., 1987. Relation of olfactory EEG to behavior: spatial analysis. Behavioral Neuroscience 101, 393–408. Freeman, W., Rogers, L., Holmes, M., Silbergeld, D., 2000. Spatial spectral analysis of human electrocorticograms including the alpha and gamma bands. Journal of neuroscience methods 95, 111–121. Galka, A., Ozaki, T., Muhle, H., Stephani, U., Siniatchkin, M., 2008. A data-driven model of the generation of human EEG based on a spatially distributed stochastic wave equation. Cognitive Neurodynamics 2, 101–113. Haykin, S., 2001. Kalman filtering and neural networks. Wiley Online Library. Hesse, W., Möller, E., Arnold, M., Schack, B., 2003. The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies. Journal of Neuroscience Methods 124, 27–44. Issa, N., Trepel, C., Stryker, M., 2000. Spatial frequency maps in cat visual cortex. Journal of Neuroscience 20, 8504. Jackson, J., 1999. Classical Electrodynamics, 3rd cd. New York, John Willey & Sons. Jansen, B., Rit, V., 1995. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics 73, 357–366. Jirsa, V., 2009. Neural field dynamics with local and global connectivity and time delay. Philosophical Transactions A 367, 1131. Jirsa, V., Jantzen, K., Fuchs, A., Kelso, J., 2002. Spatiotemporal forward solution of the EEG and MEG using network modeling. Medical Imaging, IEEE Transactions on 21, 493–504. Julier, S., Uhlmann, J., 1997. A new extension of the Kalman filter to nonlinear systems, in: Int. Symp. Aerospace/Defense Sensing, Simul. and Controls, Orlando, FL. Unscented kalman filter. 39 Julier, S., Uhlmann, J., Durrant-Whyte, H., 2002. A new approach for filtering nonlinear systems, in: American Control Conference, 1995. Proceedings of the, IEEE. pp. 1628–1632. Kaminski, M., Blinowska, K., 1991. A new method of the description of the information flow in the brain structures. Biological Cybernetics 65, 203–210. Liley, D., Alexander, D., Wright, J., Aldous, M., 1999. Alpha rhythm emerges from large-scale networks of realistically coupled multicompartmental model cortical neurons. Network: Computation in Neural Systems 10, 79–92. Ljung, L., 1999. System Identification–Theory for the user 2nd edition. Upper Saddle River, New Jersey: Prentice Hall. Lopes Da Silva, F., Blanes, W., Kalitzin, S., Parra, J., Suffczynski, P., Velis, D., 2003. Epilepsies as dynamical diseases of brain systems: Basic models of the transition between normal and epileptic activity. Epilepsia 44, 72. Lopes Da Silva, F., Van Rotterdam, A., Barts, P., Van Heusden, E., Burr, W., 1976. Models of neuronal populations: the basic mechanisms of rhythmicity. Progress in brain research 45, 281–308. Marreiros, A., Daunizeau, J., Kiebel, S., Friston, K., 2008. Population dynamics: variance and the sigmoid activation function. Neuroimage 42, 147–157. Marten, F., Rodrigues, S., Suffczynski, P., Richardson, M., Terry, J., 2009. Derivation and analysis of an ordinary differential equation mean-field model for studying clinically recorded epilepsy dynamics. Physical Review E 79, 21911. van der Merwe, R., 2003. Sigma-point kalman filters for probabilistic inference in dynamic state-space models . Moran, R., Stephan, K., Kiebel, S., Rombach, N., O’Connor, W., Murphy, K., Reilly, R., Friston, K., 2008. Bayesian estimation of synaptic physiology from the spectral responses of neural masses. NeuroImage 42, 272–284. Nunez, P., 2000. Toward a quantitative description of large-scale neocortical dynamic function and EEG. Behavioral and Brain Sciences 23, 371–398. Ozaki, T., 1993. A local linearization approach to nonlinear filtering. International Journal of Control 57, 75–96. Ozaki, T., 1994. The local linearization filter with application to nonlinear system identification , 217–240. Ozaki, T., Jimenez, J., Haggan-Ozaki, V., 2000. The role of the likelihood function in the estimation of chaos models. Journal of Time Series Analysis 21, 363–387. Peterson, D., Middleton, D., 1962. Sampling and reconstruction of wave-number-limited functions in n-dimensional Euclidean spaces. Information and Control 5, 279–323. Qubbaj, M., Jirsa, V., 2007. Neural field dynamics with heterogeneous connection topology. Physical Review Letters 98, 238102. Rasmussen, C., Williams, C., 2005. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, Cambridge MA. Rennie, C., Wright, J., Robinson, P., 2000. Mechanisms of cortical electrical activity and emergence of gamma rhythm. Journal of Theoretical Biology 205, 17–35. Robinson, P., Rennie, C., Rowe, D., O’Connor, S., 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Human brain mapping 23, 53–72. Sameshima, K., Baccalá, L., 1999. Using partial directed coherence to describe neuronal ensemble interactions. Journal of Neuroscience Methods 94, 93–103. Sanner, R., Slotine, J., 1992. Gaussian networks for direct adaptive control. IEEE Transactions on Neural Networks 3, 837–863. Sarkka, S., 2010. Continuous-time and continuous-discrete-time unscented rauch-tung-striebel smoothers. Signal Processing 90, 225 – 235. Scerri, K., Dewar, M., Kadirkamanathan, V., 2009. Estimation and model selection for an IDE-based spatio-temporal model. IEEE Transactions on Signal Processing 57, 482–492. Schiff, S., Sauer, T., 2008. Kalman filter control of a model of spatiotemporal cortical dynamics. Journal Neural Engineering 5, 1–8. 40 Schiff, S.J., 2009. Kalman meets neuron: The emerging intersection of control theory with neuroscience, in: Engineering in Medicine and Biology Society, 2009. EMBS 2009. 31st Annual International Conference of the IEEE, pp. 3318–3321. Suffczynski, P., Kalitzin, S., Lopes Da Silva, F., 2004. Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience 126, 467–484. Sullivan, T., De Sa, V., 2006. A model of surround suppression through cortical feedback. Neural Networks 19, 564–572. Ullah, G., Schiff, S., 2010. Assimilating Seizure Dynamics. PLoS Computational Biology 6, 35–45. Valdes, P., Jimenez, J., Riera, J., Biscay, R., Ozaki, T., 1999. Nonlinear EEG analysis based on a neural mass model. Biological cybernetics 81, 415–424. Wendling, F., Hernandez, A., Bellanger, J., Chauvel, P., Bartolomei, F., 2005. Interictal to ictal transition in human temporal lobe epilepsy: Insights from a computational model of intracerebral EEG. Journal of Clinical Neurophysiology 22, 343. Wilson, H., Cowan, J., 1973. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Biological Cybernetics 13, 55–80. Zetterberg, L., Kristiansson, L., Mossberg, K., 1978. Performance of a model for a local neuron population. Biological cybernetics 31, 15–26. 41 View publication stats