Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2018.00046 Applied Mathematics and Statistics Original Research Extracranial Estimation of Neural Mass Model Parameters Using the Unscented Kalman Filter Escuain-Poole Lara 1 * Garcia-Ojalvo Jordi 2 Pons Antonio J. 1 1Physics Department, Polytechnic University of Catalonia, Terrassa, Spain 2Department of Experimental and Health Sciences, Pompeu Fabra University, Barcelona, Spain

Edited by: Axel Hutt, German Meteorological Service, Germany

Reviewed by: Nicola Politi, Politecnico di Torino, Italy; Lili Lei, Nanjing University, China; Sergiy Zhuk, IBM Research (Ireland), Ireland

 *Correspondence: Lara Escuain-Poole lara.escuain@upc.edu

This article was submitted to Dynamical Systems, a section of the journal Frontiers in Applied Mathematics and Statistics

 15 10 2018 2018 4 46 22 05 2018 14 09 2018 Copyright © 2018 Escuain-Poole, Garcia-Ojalvo and Pons. 2018 Escuain-Poole, Garcia-Ojalvo and Pons

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Data assimilation, defined as the fusion of data with preexisting knowledge, is particularly suited to elucidating underlying phenomena from noisy/insufficient observations. Although this approach has been widely used in diverse fields, only recently have efforts been directed to problems in neuroscience, using mainly intracranial data and thus limiting its applicability to invasive measurements involving electrode implants. Here we intend to apply data assimilation to non-invasive electroencephalography (EEG) measurements to infer brain states and their characteristics. For this purpose, we use Kalman filtering to combine synthetic EEG data with a coupled neural-mass model together with Ary's model of the head, which projects intracranial signals onto the scalp. Our results show that using several extracranial electrodes allows to successfully estimate the state and a specific parameter of the model, whereas one single electrode provides only a very partial and insufficient view of the system. The superiority of using multiple extracranial electrodes over using only one, be it intra- or extra-cranial, is shown in different dynamical behaviours. Our results show potential toward future clinical applications of the method.

 Unscented Kalman filter data assimilation EEG neural mass model parameter estimation FI-DGR 2014-2017 2014SGR0947 FIS2015-66503 MDM-2014-0370 Generalitat de Catalunya10.13039/501100002809 Ministerio de Economía, Industria y Competitividad, Gobierno de España10.13039/501100010198 Institució Catalana de Recerca i Estudis Avançats10.13039/501100003741 1. Introduction

After several decades studying its morphology and dynamics [1], the basic mechanisms that describe the functioning of the brain are still far from being completely understood. There are different reasons that explain this arduous route toward understanding this organ. First, the neurons that form the brain are very diverse morphologically [2] and dynamically [3]. Second, these neurons are connected to each other in extremely large numbers and forming very complex networks [4], whose structural characteristics are still mostly unknown. And third, brain dynamics are very irregular and complex [5, 6]. The opposed views of an essentially noisy brain and a deterministic brain exhibiting chaotic activity have been often contrasted. On the one hand there is multiple evidence, both theoretical and experimental, that justifies a stochastic view of the brain [7, 8]. On the other hand, other studies reveal deterministic, or rather reproducible, dynamical behaviour [9, 10] both at the microscopic scale [11] and at the mesoscale recorded by electroencephalograms (EEG) or magnetoencephalograms (MEG) [12]. The reality is probably a combination of the two views. The fact that the brain receives continuous external inputs from the sensory system also makes its dynamical and experimental interpretation more complex because, even though experiments are designed to minimise uncontrolled inputs, they cannot completely rule them out. Another important limitation for studying the brain is that experimental recordings (such as EEG or fRMI) are almost always indirect reflections of the underlying neural activity [13].

A way of facing the complexities described above is by systematically comparing the experimental observations of brain activity with mathematical models based on specific hypotheses, which can thereby be validated or disproven. Modelling cerebral activity has been attempted both with top-down and bottom-up approaches [14–18]. Many of these theoretical models are simplifications that capture the basic ingredients of brain dynamics, while others are detailed accounts of the dynamics of neurons that necessarily forgo the description of the whole brain. In that context, a more feasible scale of study is the mesoscopic scale [19–25]. Many of the modern experimental techniques record information coming from populations of neurons working together. Neural mass models describe the activity of these populations mathematically using reasonably simple equations [26, 27]. These models can describe both the intrinsic oscillatory behaviour recorded at the mesoscale or event-related responses [28, 29] with morphologically plausible assumptions for their construction.

In all modelling strategies, however, identifying realistic values for the parameters of the model is a challenging task. One way to address this problem is by integrating experimental information into the models using Bayesian inference [30–35]. This strategy has started to be pursued by using Kalman filtering to integrate experimental data at both the microscopic scale of neuronal networks [36–38] and the mesoscopic scale of neural mass models [39–42]. This data assimilation approach aims to tackle the high level of noise in neuronal activity, and allows to estimate both the state and the parameters of the theoretical model using the experimental data available. The method has been used to estimate, for example, the effective connectivity that characterises epileptic seizures on a patient-specific basis (see [43] and references therein). Kalman filtering has also been used to analyse the suppression of epileptic seizures in coupled neural mass models [40, 44], and the induction of the anesthetized state by drugs [45]. But these studies use mainly invasive intracranial signals, and it would be desirable to extend them to non-invasive extracranial measurements such as EEG. Intracranial signals can be translated into EEG signals in a forward manner [46, 47], and, in the opposite direction, solving the inverse problem allows to infer intracranial signals from EEG recordings [48–50]. In this paper we advance the applications of Kalman filtering in neuroscience by extending the current procedures with a model of the head, exploring the possibilities of using non-invasive scalp measurements.

 2. Methods

To obtain a reliable estimation of the state and the dynamics of the brain, we require a biologically inspired mathematical model of its dynamics, experimental data (as non-invasive as possible), and the means of fusing both sources of information together. In this paper, for the purpose of providing a proof-of-concept of our proposed data assimilation approach, we use in silico data, instead of real experimental observations, generated by Jansen and Rit's model [26, 51], as a way to represent the dynamical evolution of the cortical structures. We then use the unscented Kalman filter as our data assimilation algorithm to estimate the state and a specific parameter of the model jointly [52–54].

2.1. Mesoscopic neural mass model

Jansen and Rit's model [26, 51] describes the mesoscopic activity of a population of neurons [55, 56], providing a good compromise between physiological realism and computational simplicity. This model reduces the neuronal diversity of a cortical column to three interacting populations: pyramidal neurons, excitatory interneurons, and inhibitory interneurons. The larger pyramidal population excites both groups of interneurons, which in turn feed back into the pyramidal cells. In our approximation, the pyramidal population is also driven by neighbouring columns and by excitatory noise representing the input from distant areas of the brain. The model is given by the following set of coupled second-order differential equations [26, 57]:

ẍ0i(t)+2aẋ0i(t)+a2x0i(t)=Aa Sigm[x1i(t)-x2i(t)], x¨1i(t)+2ax˙1i(t)+a2x1i(t)=Aa (pi(t)+k∑j=1NdKij Sigm(x1j(t−τij)                                                               −x2j(t−τij))+C2 Sigm[C1x0i(t)]), ẍ2i(t)+2bẋ2i(t)+b2x2i(t)=Bb(C4 Sigm[C3x0i(t)]),

where x0 is the average excitatory postsynaptic potential (PSP) coming to the two interneuron populations, and x1 (x2) is the average excitatory (inhibitory) PSP which inputs to the pyramidal population. The superindex i = 1 ⋯ Nd runs over all the coupled cortical columns (dipole sources) of the model. The quantity x1 − x2 is the net PSP of the pyramidal neurons, which produces the signal detected by extracranial electrodes, and is therefore our observable. The sigmoid function Sigm(v) converts the net average PSP of a population, v, into an average firing rate:

Sigm(v)=2e01+eγ(v0-v),

where e0 is the maximum firing rate of the population, γ controls the slope of the sigmoid, and v0 is the post-synaptic potential for which a 50% firing rate is obtained. The resulting firing rate is then transformed back into an average PSP by the second-order differential Equations 1–3.

The parameters A and B in the right-hand side of Equations 1–3 are the amplitudes of the excitatory and inhibitory post-synaptic potentials, and a and b are the lumped representations of the sums of the reciprocal of the time constant of the passive membrane, and all other spatially distributed delays in the dendritic network. The parameters C1 to C4 are connectivity constants that govern the interactions between populations, pi(t) is a stochastic external input that adds dynamic noise to the system, and the summation term represents the input from other coupled cortical columns. The strength of the coupling is modulated by k, with K denoting the adjacency matrix. When generating the in silico data we consider that column i receives the signal of column j with a delay τij [58]. This is because we want to generate data, in a controlled way, with a model as complex and rich in dynamics as possible to mimic real data. However, Kalman filtering, as defined in Equations 10–11 below, does not include temporal delays. Therefore, for simplicity, the model used during the filtering does not consider delays. Table 1 provides the descriptions and values of these parameters. The electrical activity detected by the electrodes on the scalp is originated by the weighted sum of the averaged membrane potential of the pyramidal cells of all the cortical columns, xi(t)=x1i(t)-x2i(t) [59], using a head model as described below.

Description and default values of the parameters for the system of neural masses.

Param. Description Value
A EPSP amplitude 3.25 mV
B IPSP amplitude 22.00 mV
a Rate constant for the excitatory population* 100 s−1
b Rate constant for the inhibitory population* 50 s−1
C1 Strength of synaptic connections from PC to EI 135
C2 Strength of synaptic connections from II to PC 108
C3 Strength of synaptic connections from PC to II 33.75
C4 Strength of synaptic connections from EI to PC 33.75
e0 Maximum firing rate of the population 2.5 s−1
v0 Mean threshold of the population 6 mV
γ Steepness of the sigmoidal transformation 0.56 mV−1
k Coupling constant 5 or 10
according to
the experiment
K Adjacency matrix Kij = 1, i ≠ j
Ki,j = 0, i = j
τ Delay Between 15
and 21 ms,
according to
distance [58]
p0 External input 200 s−1
ϵ Intensity of the noise 100 s−1

*Lumped representation of the sum of the reciprocal of the time constant of passive membrane and all other spatially distributed delays. See the section for details of the configuration of each numerical experiment. Here, PC refers to pyramidal cells, EI to excitatory interneurons, II to inhibitory interneurons, EPSP to excitatory post-synaptic potential, and IPSP to inhibitory post-synaptic potential.

 2.2. Head model

The main contribution of this paper is the use of multichannel extracranial data to obtain information about the neuronal populations inside the brain using data assimilation. To accomplish this, we use synthetic EEG data generated in silico using Jansen and Rit's model and Ary's head model. To that end, we transform the output x(t) of the neural masses to EEG signals z(t) in the electrodes (see Figure 1). This transformation is mediated by a lead field matrix [47], which builds on the basic idea of calculating the electric potential caused by a dipole source [13] on a three-layer isotropic hemisphere of radius 1 [46, 60] that represents the three main tissues that impact brain activity readings (brain, skull, and scalp). The lead field matrix also contains information about the geometry of the problem (e.g., locations of cortical columns and electrodes) and about the electrophysiology of the head (e.g., conductivities of the different tissues). The following equations show the potential Ve,i on an electrode e, located at ree [61], caused by the dipole qi(t)=xi(t)q^i generated by the cortical column i, located at rqi and oriented as q^i. In these equations, e = 1, …, Ne, where Ne is the total number of electrodes, and i = 1, …, Nd, where Nd is the total number of dipoles. Vectors are typeset in bold and modules are in regular type.

Ve,i(ree;rqi,qi)≊v1(ree;μ1rqi,ρ1qi)+v2(ree;μ2rqi,ρ2qi)+v3(ree;μ3rqi,ρ3qi), v1(ree;rqi,qi)=((c1e,i,1−c2e,i,1(ree·rqi))rqi+                                                                                    c2e,i,1(rqi)2ree)·qi, v2(ree;rqi,qi)=((c1e,i,2−c2e,i,2(ree·rqi))rqi                                                                                  +c2e,i,2(rqi)2ree)·qi, v3(ree;rqi,qi)=((c1e,i,3−c2e,i,3(ree·rqi))rqi                                                                                  +c2e,i,3(rqi)2ree)·qi.

In these expressions,

c1e,i,s=14πσs(rqi)2(2de,i·rq i(de,i)3+1de,i−1ree),c2e,i,s=14πσs(rqi)2(2(de,i)3+de,i+reereΓ(re e,rq i)),Γ(re e,rq i)=de,i(reede,i+(ree)2−(rq i·re e)).

Extracranial data generation and illustration of Ary's model of the head. The light and dark red arrows indicate dipole sources, and the electrodes are shown as gray and black rectangles. The elements in the cartoon illustrate how all the signals produced by the cortical columns (represented with the solid red line in the top left panel) are transformed into an electrode reading (shown in black dots in the top right panel) through the lead field matrix. In this drawing, as in Equations (5–9), rq is the distance from the origin to the cortical column under consideration; re is the distance from the origin to the electrode; and d is the distance from the cortical column to the electrode. The placement of the arrows here is for illustration purposes only; in our study, the cortical columns are placed on the surface of the brain, close to the skull.

The tangential conductivity of each layer is represented by σs [60] and ρs and μs are the Berg parameters relative to it [62] (see Table 2). The parameter de,i=re e−rq i is the relative position of the electrode e under consideration with respect to the position of the dipole i.

Values of the Berg parameters for the three layers [60, 62].

Parameter Layer 1 Layer 2 Layer 3
Tangential conductivity σs 1.0 0.0125 1.0
Berg parameter ρs 0.9901 0.7687 0.4421
Berg parameter μs 0.0659 0.2389 0.3561

Cartesian coordinates of the dipoles used throughout the study.

x y z
dipole 1 0.1688 0.2242 0.2597
dipole 2 0.3766 –0.8520 0.2597
dipole 3 0.6622 –0.2242 –0.1948

The origin of coordinates is the centre of the perimeter of the head.

 2.3. The unscented Kalman filter for data assimilation

The Unscented Kalman Filter (UKF) is our algorithm of choice to bring together the dynamical state of the model and the in silico data. It is a standard tool in the field of systems and control engineering, and has been shown to be both computationally efficient and robust even when dealing with stochastic nonlinear systems [63]. In our case, the computational burden—O(n3), where n is the size of the state—is acceptable for a biologically reasonable number of sources. In order to simultaneously estimate the state and parameters of the model described by Equations (1)–(3), we regard it as a discrete-time state-space dynamical system of the following form:

xk+1=F(xk)+vk zk=H(xk)+wk

where x=(x01,x11,x21,x02,…,x2Nd,θ1,…,θNp)∈ℝnx is the state vector (related to the variables and parameters of the model), with θp being the parameters to estimate, which obey the equations θ˙p=0. (In our joint estimation of the parameters, these are included in the state vector together with the system variables). The vector z∈ℝnz is the measurement vector (our in silico EEG readings). The vectors v and w are uncertainty terms that account for process noise and measurement noise, respectively, with Gaussian distributions p(v) ~N(0, Q) and p(w) ~N(0, R), respectively. The process transition F is obtained with a numerical implementation of Equations (1)–(3), as described below. Finally, H relates the state to measurement space, which is either Ve, i (in the case of simulated EEG), or xi (in the case of simulated electrocorticography). Interestingly, where EEG is concerned, this basic part of the Kalman filter is in our case implemented by the skull, the effect of which is represented by the lead field matrix, based on Ary's head model and introduced above.

The UKF is a recursive predictor-corrector-type algorithm that aims to minimise the mean square error of the estimated states and parameters over time. For each time step it calculates a prediction of the state and parameters of the system, which is corrected when the information from a measurement is incorporated. The amount of confidence given to the model and measurement is quantified by the Kalman gain K, which is calculated at each time step based on prediction covariances as well as model and measurement error covariances (Q and R, respectively). For more details on the implementation of the filter, the reader is referred to the Appendix and to Kalman [54], Merwe and Wan [52], Julier and Uhlmann [53], and Solonen et al. [64].

 2.4. Generation of in silico datasets

For this paper three different in silico datasets were generated. We consider both simulated electrocorticography (ECoG, intracortical) and electroencephalography (EEG, extracranial) readings (using Ary's model in the latter case). We chose to use three sources because this provides a considerable spatial and temporal richness in the resulting signals, while keeping the system reasonably simple and still biologically plausible [65, 66].

The presence of additional dipoles in the brain, and its influence on the sources of study, is accounted for in the stochastic external input to the sources (p(t), see Equation 2):

p(t)=p0+ξ(t),

where p0=200 s-1 and ξ(t) is Gaussian white noise [67] of zero mean and correlation 〈ξ(t)ξ(t′)〉 = 2ϵδ(t − t′) [68]. At the extracranial level, the other sources also affect the final EEG signal, as well as the different tissues (brain, skull, scalp, and even hair). This is modelled by adding Gaussian noise with zero mean and standard deviation 100 mV (unless otherwise stated) to the simulated EEG.

All datasets used the same locations for the cortical columns [66]. The electrodes were placed using a subset of the equidistant layout, a standard layout for EEG [69] (roughly illustrated in Figures 6–8). The strength of the coupling was set at a medium value so that the cortical columns have a visible effect on one another without fully synchronising behaviours and locking their dynamics (between k = 5 and k = 10), and the configurations of the couplings are as shown in Figure 2. Table 1 shows representative values for the parameters used in all analyses unless otherwise specified. In this paper we focus on estimating the amplitudes A of the EPSPs of the different cortical columns, and therefore we choose values for these amplitudes that produce signals that reflect various dynamic regimes that we wish to explore. (The rest of the parameters were fixed to their standard values [26, 51], as described in Table 1).

The two cortical column motifs used in this paper. Unidirectionally coupled cortical columns have no backflow (A), and bidirectionally coupled columns are coupled all-to-all (B). See Table 3.

The numerical solver used to generate the in silico time series was the Heun algorithm [70] with a time step of Δt = 1 ms. The length of the data is 100 s in all cases. Using the Heun algorithm together with Equations 1– 3 to update the state variables and the lead field matrix (in order to get the potential in the electrodes of the scalp in Equations 5– 9), we generate the required map to apply Kalman filtering in Equations 10 and 11. The following equations implement the stochastic Heun algorithm used to update xk:

xk+1=xk+12(F(xk)+F(x~k))Δt+12∑(g(xk)+g(x~k))X, x~k=xk+F(xk)Δt+g(xk)X.

Where g(…), together with Equation 12, introduces the noise term in Equation 2 and is zero for Equations 1 and 3. In X=2ϵΔtγ, γ are gaussianly distributed random numbers with zero mean and unit variance. At different instants of time, these random numbers are independent from one another.

2.4.1. Three unidirectionally coupled cortical columns

For the first study the cortical columns were coupled unidirectionally (Figure 2A), as described in Liu and Gao [71]. The parameters were set to standard values [26] for the three cortical columns (see Table 1), except for the first column, in which A1 was set to 3.58 mV to make it hyperexcitable. Additionally, the three cortical columns had p0 = 90 s−1 and ϵ = 2 s−1. This first hyperexcitable column causes a spiking cascade in the other two columns. With this experiment, we aimed to compare how extra- and intra-cranial electrodes perform in the case of a behaviour being induced by an input from another column, and not by the column's own parameter configuration. The resulting data can be found in the Data Sheets 1, 2 in Supplementary Material. Please, see the README file (Data Sheet 6) for more information.

 2.4.2. Three bidirectionally coupled cortical columns: coarse parameter estimation

The three cortical columns are located as in the previous section, but coupled bidirectionally (Figure 2B). Additionally, the maximum amplitudes of the excitatory PSPs were set to A1 = 4.25 mV, A2 = 10.00 mV, and A3 = 3.25 mV. These values were chosen to cause the three cortical columns to be in very different dynamical regimes: cortical column 1 operates in a spiking regime; cortical column 2 oscillates with alpha frequency but with an amplitude similar to that of the spikes; and cortical column 3 oscillates in a more standard regime, as described in [26]. Also, the external input p(t) for each of the three cortical columns was set using p0=200 s-1 and ϵ = 100 s−1. Our aim here was to study how the filter performs in an extreme situation, in which the dynamics of the columns are widely different from one another. We intended to explore the outcome of estimating with single extracranial electrodes as well as the complete set, and to compare with intracranial estimation (Data Sheets 3, 4).

 2.4.3. Three bidirectionally coupled cortical columns: fine parameter estimation

In the previous section, the value of A of one of the cortical columns was much larger than the other two. We now consider the same coupling motif, but with values of the A parameter that are much closer together in value: A1 = 3.58 mV, A2 = 3.25 mV, and A3 = 3.10 mV. (The values defining the external input p(t) remain the same as in the previous experiment). Our goal was to check if the filter can discriminate between the values when they are closer together (Data Sheet 5).

 2.5. Filtering

For each of the experiments we conducted 50 realisations of each estimation for the complete state vector, with different initial conditions; all the figures show averages of the 50 estimations, unless otherwise specified. The initial conditions for state and parameter estimations were randomly generated with a normal distribution of zero mean and unit variance; the parameters, however, were constrained to deviate no more than 90% of their actual value as an initial assumption.

The noise covariances Q and R were chosen according to the best knowledge of the system and of the noise corrupting the data. Therefore, Q was set to contemplate the incoming noise to each dipole, i.e., it was set to a null matrix except for the term corresponding to the equation that contains the input p(t) (see Equation 2 and [68]). The matrix R was set to 1000I mV2. (In practice, in most applications of the Kalman filter, the matrix R is fairly easy to set with the knowledge of the measurement precision as a starting point, but Q is often set by trial and error).

 2.6. Ethics statement

All data used in this manuscript come from numerical simulations of a mathematical model. No human or animal data have therefore been used, and ethics approval was not necessary.

 3. Results

In order to compare the performance of the extra- and intra-cranial approaches to Kalman filtering, we have analysed three different cortical column configurations, each using one of the two motifs shown in Figure 2. Where relevant, two different types of estimations have been used: intracranial and extracranial. Intracranial estimation uses simulated data that would have hypothetically been obtained from electrocorticography, that is, using a single intracortical electrode, and is estimated with the data provided by a single location—in other words, the direct output of Jansen and Rit's model. Extracranial estimation, on the other hand, employs simulated data originated from EEG recordings, using several electrodes placed on the skull, and is implemented here with the projection on the head of the model output. We now discuss the results for the three different situations that we have considered.

3.1. Three unidirectionally coupled cortical columns

In this case, information flows unidirectionally because of the way the cortical columns are coupled [71]. As can be seen in the lower panels of Figure 3, the first cortical column has a random spiking activity, due to the increased value of A and the presence of noise [20]. Due to the architecture of the coupling, cortical column 1 causes cortical columns 2 and 3 to spike also, when otherwise they would have simply fluctuated around their resting level.

Intracranial and extracranial fittings with propagated excitation along unidirectionally coupled cortical columns. The upper panels show the estimation of parameter A, and the lower panels show the estimations of the observed states. The lines show the averages of the 50 realisations of the estimation, and the shadowed areas indicate the standard deviation. The actual values of the A parameters are A1 = 3.58 mV, A2 = 3.25 mV, and A3 = 3.25 mV, the other parameters being set to standard values (Table 1). All three cortical columns received an external input, p(t), with p0=90 s-1 and ϵ = 2 s−1. The coupling constant was set to k = 10. The measurements were corrupted with Gaussian noise of mean 0 and standard deviation 100 mV for extracranial measurements and standard deviation 5 mV for intracortical measurements. Except for cortical column 1, with intracortical data the filter converges to a much higher value than the target, whereas with extracranial data the filter converges to a value which is accurate. In the lower panels it is shown that extracranial estimations of the state are also accurate, whereas intracortical estimations fail to reproduce the spikes correctly.

The upper panels of Figure 3 show the intracortical and extracranial estimations of A for the three cortical columns. The estimation for A1 of the first column converges to its correct value, with both the intra- and extracortical approaches. This was to be expected, since the first cortical column receives no inputs from other elements of the system. In contrast, the intracortical estimations for cortical columns 2 and 3 converge to values significantly higher than their actual value of 3.25 mV. We conjecture that this is caused by the spiking of these two cortical columns, which as mentioned above is due to the influence of cortical column 1. Multi-channel extracranial information, in contrast, allows to see the complete picture of the coupled cortical columns and treat them as a single composed system, contrary to the partial picture obtained from the information provided by the single intracranial recordings. Therefore, estimation is better when using extracranial information with several electrodes, as shown in the upper panels of the figure. The lower panels of Figure 3 show the estimation of the state. The UKF shows great efficacy when the estimation is extracranial, but performs poorly in the case of intracortical estimation (with the exception of cortical column 1, because it has no input from other cortical columns). This highlights the value of extracranial estimation, in which it is possible to take the whole brain into account in a non-invasive manner.

 3.2. Three bidirectionally coupled cortical columns: coarse parameter estimation

The second experiment aims to explore the possibilities of the filter in more extreme situations, as the parameters were chosen to reflect more diverse dynamical regimes. The following sections describe the results of single- and multichannel estimations.

3.2.1. Moderate intracortical measurement noise

Figure 4 shows again the performance obtained using the simulated data from a set of extracranial electrodes compared to using individual intracortical electrodes for each cortical column. In this case we show the 50 realisations of each filtering, without showing the average. The extracranial data for this experiment were corrupted with a measurement Gaussian noise of zero mean and standard deviation 100 mV; the intracortical data were corrupted with a measurement noise of standard deviation 5 mV in order to maintain similar levels of signal-to-noise ratio.

Intracranial and extracranial fittings for coarse parameter estimation in the case of bidirectional coupling. As in the previous figure, the upper panels show the estimation of A for each cortical column and the lower show the estimations of the observed states. The results are shown here without averaging. The actual values of the amplitudes of the EPSPs are A1 = 4.25 mV, A2 = 10.00 mV, and A3 = 3.25 mV; the rest of the parameters were set to standard values (Table 1). The external input p(t) for the three cortical columns had p0 = 200 s−1 and ϵ = 100 s−1. The coupling constant was set to k = 5. The intracortical measurements were corrupted with Gaussian noise of mean 0 and standard deviation 5 mV, while the noise in the extracranial measurements has standard deviation 100 mV. Extracranial estimations of the parameters are both faster and more accurate than intracortical estimations; this applies also to the state, whose dynamics are more faithfully reproduced using multi-electrode extracranial estimation (as shown in the lower panels).

As shown in Figure 4, the intracortical parameter estimations do not approximate the target value very well. In particular, the estimations of A for cortical column 2 converge to three different values depending on the initial conditions. The state estimation follows the actual state of the system closely only for cortical column 1. The situation is very different when with extracranial electrodes, where all 50 realisations of the estimations converge with much more precision to the correct values for both state and parameters (with the exception of A2, which still tends to lower values in a very small quantity of the realisations). Again, extracranial performance is better, in general, to intracortical.

 3.2.2. High intracortical measurement noise

The difference between intracranial and extracranial estimation is even larger for higher measurement noise (Figure 5). In this case, the amount of noise in the intracortical data was set to the same value as the noise in the extracranial data. The value of R was tuned to reflect the increase in measurement noise, but the intracortical estimations failed to obtain the correct values for the parameters and reproduce the state.

Intracranial and extracranial fittings for coarse parameter estimation, with a higher amount of intracortical measurement noise. The upper panels show the estimation of the EPSPs for each cortical column and the lower panels show the estimations of the observed states. The results are shown here without averaging. The actual values of the amplitudes of the EPSPs are A1 = 4.25 mV, A2 = 10.00 mV, and A3 = 3.25 mV; the rest of the parameters were set to standard values (Table 1). The external input p(t) for the three cortical columns had p0=200 s-1 and ϵ = 100 s−1. The coupling constant was set to k = 5. The intracortical measurements were corrupted with Gaussian noise of mean 0 and standard deviation 100 mV—about an order of magnitude higher than the noise in the previous graph—, while the noise in the extracranial measurements has standard deviation 100 mV. Extracranial estimations of the parameters are also faster and more accurate than intracortical estimations, more markedly so in this case; as to the state, in this more extreme case, the intracortical estimation does not mimic the evolution of the system in any way.

 3.2.3. Using one single extracranial electrode

Using the same dataset, we aimed to investigate the outcome of using each extracranial electrode individually [43], as opposed to using the complete subset as until now. Therefore, we used each electrode separately to estimate the state and parameters of the complete system, with 50 realisations of the estimation for each electrode. By doing so, we show that the quality of the estimations is strongly dependent on the relative positions of sources and electrodes.

In Figures 6–8 we present the results for the estimation of parameter A of each of the three cortical columns separately. The histograms show the distribution of the 50 estimations of A using each electrode, placed in the respective position of the electrode in question. Vertical coloured lines in the histograms mark the value of the three A parameters being estimated (one in each figure). The histograms show a strong dependence on space of the quality of the estimations. As a general trait, the estimations are better when the electrodes are near the cortical column whose value of A is being estimated, whereas the more distant electrodes show a wider distribution of final values for the parameter.

Distribution of 50 realisations of A estimations from a single electrode for cortical column 1 (solid red circle). The histograms are placed at the location of the corresponding measuring electrode, and the location of the three cortical columns generating the activity are shown with coloured circles (with the full circle corresponding to the column whose value of A is being estimated in this figure). Vertical lines with the same colours as the circles mark the corresponding actual A values. The distributions tend to be narrowest in the vicinities of cortical column 1. Nevertheless, they do not group around the target value of A1 = 4.25 mV (vertical red line), as they should, but around that of A3 = 3.25 mV (vertical blue line).

In Figure 6 the distribution of the estimations of A1 are shown. The distributions tend to be narrowest in the vicinities of the cortical column whose A value is being estimated. However, it is noteworthy that the histograms obtained from the observations in distant electrodes tend to group not around the actual value of A1 = 4.25 mV (red vertical line), but of A3 = 3.25 mV (blue vertical line). This result suggests that the algorithm is unable to distinguish the origin of the EEG activity when sources and electrodes are distant from each other.

Figure 7 shows the results of the estimation of A2 (actual value shown by vertical green lines), revealing wider distributions in general, which indicates a stronger dependence on initial conditions. Although it is true that the electrodes near cortical column 2 perform better in estimating A for that column, the difference with more distant electrodes is not as large as for the estimates of A for cortical columns 1 and 3.

Distribution of 50 realisations of A estimations from a single electrode for cortical column 2 (solid green circle). The distributions here are wider than for A1 and A3, although they still tend to be more accurate near the cortical column (solid green circle) and group around the target value of A2 = 10.00 mV (vertical green line).

Finally, Figure 8 shows the performance of each electrode when A3 is being estimated (actual value shown by vertical blue lines in the figure). Interestingly, even the electrodes located at the far left of the figure lead to a good estimate of A, comparable to that coming from the electrodes in the far right, which are closer to column 3 and could therefore be expected to provide a much more accurate estimation.

Distribution of 50 realisations of A estimations from a single electrode for cortical column 3 (solid blue circle). As in the two previous figures, the distributions for the electrodes closest to the source (solid blue circle) are narrow, grouping around the correct value (A3 = 3.25 mV, vertical blue line). Surprisingly, the electrodes in the far left also give rise to narrow distributions.

While the estimations arising from single electrodes are reasonably accurate in some cases, using the complete set of 15 electrodes invariably yields better results. This is because, in Kalman filtering, combining many sources of information always improves the final estimation, even if some of the sources are inaccurate or incomplete [72].

 3.3. Three bidirectionally coupled cortical columns: fine parameter estimation

In the previous section, the aim was to generate widely different dynamics in each column. We now consider the results of estimating parameters which are much closer to one another. The purpose of this test was to ascertain whether the filter could differentiate between parameters with smaller differences in value. This ability is very important if we expect to use the technique in clinical applications. Figure 9 shows the extracranial estimation of the A parameters using the complete subset of 15 electrodes. The estimations converge to the actual values with enough accuracy as to give hopes of using the filter in a clinical setting.

Extracranial fit with parameters close together in value. The estimations of the amplitude of the EPSPs of the three cortical columns are shown after averaging over 50 realisations (solid lines); the shadowed areas indicate the standard deviation. The actual values of the amplitudes of the EPSPs are A1 = 3.58 mV, A2 = 3.25 mV, and A3 = 3.10 mV; the rest of the parameters were set to standard values (Table 1). The external input p(t) for the three cortical columns had p0=200 s-1 and ϵ = 100 s−1. The coupling constant was set to k = 5. The Gaussian noise in the extracranial measurements has standard deviation 100 mV. The estimation of the parameters is fairly accurate.

 4. Discussion

The most important limitation of current data assimilation processes in neuroscience is that the appropriate experimental recordings are usually intracranial. Despite this fact, using Kalman filtering to fit these data to neural mass models shows promise in several contexts and applications. In this study we have modified this type of approach by extending it with a head model, with the aim of integrating non-invasive experimental recordings taken from the scalp (EEG). By increasing the range of recordings in this way, the application of the data assimilation protocols opens up to the large set of situations in which scalp recordings are used. We keep the exploration of the technique using real EEG experimental data in mind, but here we have explored the limitations and advantages of our model using in silico data in very well controlled conditions.

Our main goal in this paper has been to show that data assimilation employing multiple non-invasive EEG electrodes (as coming from scalp EEG measurements) provides a better estimate of the brain's dynamical state than using a single invasive (intracranial) EEG electrode. In particular, we have aimed at contrasting our results with existing work using the latter approach, which has employed a filtering method, namely the unscented Kalman filter [39]. Filtering methods have been so far the method of choice in data assimilation problem in neuroscience [37, 38, 40–43], with variational methods having been used very sparsely [73]. We thus chose to work with a filtering algorithm, the UKF, that is already relatively well characterized in neural model, and which we could therefore use as a benchmark.

We have considered a system comprised of three cortical columns, modelled according to Jansen and Rit's equations and coupled following two different motifs. The cortical columns are all driven by a noisy input coming from the columns of the rest of the brain and sensory stimuli. The signal from the cortical columns is then transferred to the skull, after which it is corrupted with Gaussian noise to simulate electrode readings from EEG. These are then used to estimate the amplitude of the excitatory post-synaptic potentials.

Even though the quality of the experimental measurements at the scalp might be, in general, worse than the intracranial recordings, EEG can always be measured from several positions. This allows to obtain measurements for patients without intracranial implants and also to compensate the potentially low quality of the data by having many recordings at the same time. Besides, the spatial distribution of the electrodes on the scalp allows the information arriving from the whole cortex to be available during the assimilation process. In order to address these strengths and weaknesses of the scalp recordings with respect to intracranial measurements, we have analysed situations where assimilation with only intracortical recordings may be wanting, where diverse dynamical regimes coexist due to large differences in control parameters in the cortical columns, or where fine changes of the parameters make the discrimination difficult.

The first study considered here involves three columns that are coupled unidirectionally with no backflow. The first cortical column is made hyperexcitable by increasing the excitatory post-synaptic potential to A1 = 3.58 mV; this cortical column causes the second cortical column and, indirectly, the third, to modify their behaviour by inducing spiking. For the intracranial estimations, single intracortical electrodes measured the evolution of the three cortical columns independently; for the extracranial estimations, 15 extracranial electrodes were used simultaneously. Applying the Kalman filter to the extracranial data provided a good estimation of the A parameters and of the dynamical state of the model; the intracortical measurements, however, yielded mixed results. The estimation for cortical column 1 was accurate, whereas for cortical columns 2 and 3 the estimation of A was above the target value and very close to the estimation for cortical column 1 (see orange dashed lines in Figure 3). The estimation of the dynamical state of cortical columns 2 and 3 was also worse than the estimation for cortical column 1. We attribute this to the fact that columns 2 and 3 are excited by column 1, which spikes due to a higher value of A. As a consequence, when independently evaluated using the intracranial information, the estimation is higher than the actual value. Therefore we suggest that one intracranial electrode provides only a partial view of the system, and thus cannot capture the behaviours of all three cortical columns and the interactions between them; the use of many electrodes provides a more complete view of the system.

Next we considered a situation in which the dipoles were coupled bidirectionally in an all-to-all configuration. The A parameters were chosen such as to cause different dynamic behaviours in the cortical columns. Three types of fitting via Kalman filtering were performed, using (i) independent intracortical recordings of single cortical columns, (ii) the complete subset of 15 extracranial electrodes, and (iii) single extracranial electrodes. The intracortical data were corrupted with two different levels (medium and high) of measurement noise. For both cases, the multi-electrode extracranial estimation surpasses the intracortical results in both speed of convergence and quality; the difference, however, is more marked in the presence of higher measurement noise in the intracortical recordings. In all these cases, the representation of the dynamical state of the three cortical columns using the complete set of 15 extracranial electrodes nicely matched the actual dynamical state, contrary to the limited match obtained using single intracranial or extracranial recordings. The results for the single electrodes show a significant influence of space on the quality of the estimations, in the sense that estimations of electrodes close to the source are relatively accurate, and electrodes further away from the source might not allow to discriminate the source of the information correctly, or might completely fail to represent the system.

Finally, we considered the situation of an identical cortical column configuration—in terms of situation and coupling—, except for the values of the EPSPs of the cortical columns. This dataset was filtered only extracranially, with the purpose of evaluating the filter's ability to discriminate parameter values within narrower ranges. The results in this case were also reasonably good, even though the real values of the parameter were much closer to one another, which makes data assimilation more challenging.

Even though the results shown here are better when considering extracranial electrodes, the method has, of course, limitations. For instance, the head model introduces new parameters which should be realistic. The use of Jansen's model, while being a very standard choice in the field, is not mandatory and could be substituted by others. There several alternatives to Ary's head model too. The succesful application of the method with different combinations of these models will, for sure, guide researchers to choose which models are more suitable for the theoretical description of the mesoscale in the brain. Even though the exploration of the dynamics for the different neural mass models or of the different head models might be worth exploring in future works, it lays outside of the scope of this work.

Applications of the method presented here will certainly appear in the field of brain-machine interface, long-term tracking for early diagnosis of degenerative diseases, or short-term tracking during rehabilitation of traumas and strokes. However, the succesful application of the method in each of these fields will require further research.

Taken as a whole, our results show that, independently of the need to explore more realistic situations, extracranial EEG recordings constitute a good candidate to be used together with neural mass models and Kalman filters, provided the method is extended with a head model. With its management of the noise in the system and of the inherent simplifications in neurological models, the Kalman filter is an appropriate tool for tackling the challenges of brain data processing. Using non-invasive techniques in these processes widens the applications of Kalman-based data assimilation methods in neuroscience.

 Author contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors wish to thank Drs. Maciej Jedynak and Daniel Malagarriga for helpful discussions.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2018.00046/full#supplementary-material

 References Yuste R. From the neuron doctrine to neural networks. Nature Rev Neurosci. (2015) 16:487–97. 10.1038/nrn396226152865 Braitenberg V Schütz A. Anatomy of the cortex. In: Statistics and Geometry, No. 18 in Studies of Brain Function. Berlin: Springer (1991). Izhikevich EM Edelman GM. Large-scale model of mammalian thalamocortical systems. Proc Natl Acad Sci USA. (2008) 105:3593–8. 10.1073/pnas.071223110518292226 Callaway EM. Micro-, Meso- and Macro-Connectomics of the Brain. Springer (2016). Available online at: http://link.springer.com/10.1007/978-3-319-27777-6 Freeman WJ. Simulation of chaotic EEG patters with a dynamic model of the olfactory system. Biol Cybern. (1987) 56:139–50. 10.1007/BF003179883593783 Faure P Korn H. Is there chaos in the brain? I. Concepts of nonlinear dynamics and methods of investigation. C R Acad Sci III. (2001) 324:773–93. 10.1016/S0764-4469(01)01377-411558325 Shadlen MN Newsome WT. Noise, neural codes and cortical organization. Curr Opin Neurobiol. (1994) 4:569–79. 10.1016/0959-4388(94)90059-07812147 Faisal AA Selen LPJ Wolpert DM. Noise in the nervous system. Nature Rev Neurosci. (2008) 9:292–303. 10.1038/nrn225818319728 Schiff SJ Jerger K Duong DH Chang T Spano ML Ditto WL. Controlling chaos in the brain. Nature (1994) 370:615–20. 10.1038/370615a08065447 van Vreeswijk C Sompolinsky H. Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science (NY) (1996) 274:1724–6. 10.1126/science.274.5293.17248939866 Celletti A Froeschle C Tetko IV Villa AEP. Deterministic behaviour of short time series. Meccanica. (1999) 34:147–54. 10.1023/A:1004668310653 Stam CJ. Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. Clin Neurophysiol. (2005) 116:2266–301. 10.1016/j.clinph.2005.06.01116115797 Buzsáki G Anastassiou Ca Koch C. The origin of extracellular fields and currents–EEG, ECoG, LFP and spikes. Nat Rev Neurosci. (2012) 13:407–20. 10.1038/nrn324122595786 Wright JJ Liley DTJ. Dynamics of the brain at global and microscopic scales: neural Networks and the EEG. Behav Brain Sci. (1996) 19:285. 10.1017/S0140525X00042679 Rabinovich MI Varona P Selverston AI Abarbanel HDI. Dynamical principles in neuroscience. Rev Modern Phys. (2006) 78:1213. 10.1103/RevModPhys.78.1213 Eliasmith C Stewart TC Choo X Bekolay T DeWolf T Tang C . A large-scale model of the functioning brain. Science (2012) 338:1202–5. 10.1126/science.122526623197532 Deco G Ponce-Alvarez A Hagmann P Romani G Mantini D Corbetta M. How local excitation-inhibition ratio impacts the whole brain dynamics. J Neurosci. (2014) 34:7886–98. 10.1523/JNEUROSCI.5068-13.201424899711 Nogaret A Meliza CD Margoliash D Abarbanel HDI. Automatic construction of predictive neuron models through large scale assimilation of electrophysiological data. Sci Rep. (2016) 6:32749. 10.1038/srep3274927605157 David O Friston KJ. A neural mass model for MEG/EEG: coupling and neuronal dynamics. Neuroimage (2003) 20:1743–55. 10.1016/j.neuroimage.2003.07.01514642484 Grimbert F Faugeras O. Bifurcation analysis of Jansen's neural mass model. Neural Comput. (2006) 18:3052–68. 10.1162/neco.2006.18.12.305217052158 Cona F Zavaglia M Massimini M Rosanova M Ursino M. A neural mass model of interconnected regions simulates rhythm propagation observed via TMS-EEG. Neuroimage (2011) 57:1045–58. 10.1016/j.neuroimage.2011.05.00721600291 Coombes S. Large-scale neural dynamics: Simple and complex. Neuroimage (2010) 52:731–9. 10.1016/j.neuroimage.2010.01.04520096791 Babajani A Soltanian-Zadeh H. Integrated MEG/EEG and fMRI model based on neural masses. IEEE Trans Biomed Eng. (2006) 53:1794–1801. 10.1109/TBME.2006.87374816941835 Babiloni F Cincotti F Babiloni C Carducci F Mattia D Astolfi L . Estimation of the cortical functional connectivity with the multimodal integration of high-resolution EEG and fMRI data by directed transfer function. Neuroimage (2005) 24:118–31. 10.1016/j.neuroimage.2004.09.03615588603 Bojak I Oostendorp TF Reid AT Kötter R. Connecting mean field models of neural activity to EEG and fMRI data. Brain Topogr. (2010) 23:139–49. 10.1007/s10548-010-0140-320364434 Jansen BH Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. (1995) 73:357–66. 10.1007/BF001994717578475 Wendling F Bartolomei F Bellanger JJ Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Euro J Neurosci. (2002) 15:1499–508. 10.1046/j.1460-9568.2002.01985.x12028360 David O Harrison L Friston KJ. Modelling event-related responses in the brain. Neuroimage (2005) 25:756–70. 10.1016/j.neuroimage.2004.12.03015808977 Spiegler A Knösche TR Schwab K Haueisen J Atay FM. Modeling brain resonance phenomena using a neural mass model. PLoS Comput Biol. (2011) 7:e1002298. 10.1371/journal.pcbi.100229822215992 Wright JJ Robinson PA Rennie CJ Gordon E Bourke PD Chapman CL . Toward an integrated continuum model of cerebral dynamics: the cerebral rhythms, synchronous oscillation and cortical stability. BioSystems (2001) 63:71–88. 10.1016/S0303-2647(01)00148-411595331 Schiff SJ Sauer T. Kalman filter control of a model of spatiotemporal cortical dynamics. J Neural Eng. (2008) 5:1–8. 10.1088/1741-2560/5/1/00118310806 Kiebel SJ Garrido MI Moran RJ Friston KJ. Dynamic causal modelling for EEG and MEG. Cogn Neurodyn. (2008) 2:121–36. 10.1007/s11571-008-9038-019003479 Shine JM Koyejo O Bell PT Gorgolewski KJ Gilat M Poldrack RA. Estimation of dynamic functional connectivity using Multiplication of Temporal Derivatives. Neuroimage (2015) 122:399–407. 10.1016/j.neuroimage.2015.07.06426231247 Ma X Chou CA Sayama H Chaovalitwongse WA. Brain response pattern identification of fMRI data using a particle swarm optimization-based approach. Brain Informat. (2016) 3:181–92. 10.1007/s40708-016-0049-z27747594 Alswaihli J Potthast R Bojak I Saddy D Hutt A. Kernel reconstruction for delayed neural field equations. J Math Neurosci. (2018) 8:3. 10.1186/s13408-018-0058-829399710 Voss HU Timmer J Kurths J. Nonlinear dynamical system identification from uncertain and indirect measurements. Int J Bifurcat Chaos (2004) 14:1905–33. 10.1142/S0218127404010345 Hamilton F Berry T Peixoto N Sauer T. Real-time tracking of neuronal network structure using data assimilation. Phys Rev E. (2013) 88:052715. 10.1103/PhysRevE.88.05271524329304 Hamilton F Cressman J Peixoto N Sauer T. Reconstructing neural dynamics using data assimilation with multiple models. EPL (Europhys Lett). (2014) 107:68005. 10.1209/0295-5075/107/68005 Freestone DR Kuhlmann L Chong MS Nesic D Grayden DB Aram P . Patient-specific neural mass modeling - stochastic and deterministic, methods. In: Recent Advances in Predicting and Preventing Epileptic Seizures. (2013) p. 63–82. Available online at: https://hal.archives-ouvertes.fr/hal-00876475/%5Cnhttp://www.worldscientific.com/doi/abs/10.1142/9789814525350_0005 Shan B Wang J Deng B Wei X Yu H Li H. UKF-based closed loop iterative learning control of epileptiform wave in a neural mass model. Cogn Neurodyn. (2015) 9:31–40. 10.1007/s11571-014-9306-026052360 López-Cuevas A Castillo-Toledo B Medina-Ceja L Ventura-Mejía C. State and parameter estimation of a neural mass model from electrophysiological signals during the status epilepticus. Neuroimage. (2015) 113:374–86. 10.1016/j.neuroimage.2015.02.05925754070 Li S Li J Li Z. An improved unscented kalman filter based decoder for cortical brain-machine interfaces. Front Neurosci. (2016) 10:587. 10.3389/fnins.2016.0058728066170 Freestone DR Karoly PJ Nešić D Aram P Cook MJ Grayden DB. Estimation of effective connectivity via data-driven neural modeling. Front Neurosci. (2014) 8:383. 10.3389/fnins.2014.0038325506315 Cao Y Ren K Su F Deng B Wei X Wang J. Suppression of seizures based on the multi-coupled neural mass model. Chaos (2015) 25:103120. 10.1063/1.493171526520086 Kuhlmann L Freestone DR Manton JH Heyse B Vereecke HEM Lipping T . Neural mass model-based tracking of anesthetic brain states. Neuroimage (2016) 133:438–56. 10.1016/j.neuroimage.2016.03.03927018048 Zhang Z. A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres. Phys Med Biol. (1995) 40:335–49. 10.1088/0031-9155/40/3/0017732066 Mosher JC Leahy RM Lewis PS. EEG and MEG: forward solutions for inverse methods. IEEE Trans Biomed Eng. (1999) 46:245–59. 10.1109/10.74897810097460 Haufe S Tomioka R Dickhaus T Sannelli C Blankertz B Nolte G . Large-scale EEG/MEG source localization with spatial flexibility. Neuroimage (2011) 54:851–9. 10.1016/j.neuroimage.2010.09.00320832477 Verhellen E Boon P. EEG source localization of the epileptogenic focus in patients with refractory temporal lobe epilepsy, dipole modelling revisited. Acta Neurol Belg. (2007) 107:71–7. Available online at: http://www.actaneurologica.be/pdf.php?annee=2007&trimestre=3 18072334 Gotman J. Noninvasive methods for evaluating the localization and propagation of epileptic activity. Epilepsia. (2003) 44(Suppl 1):21–9. 10.1111/j.0013-9580.2003.12003.x14641558 Jansen BH Zouridakis G Brandt ME. A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol Cybern. (1993) 68:275–83. 10.1007/BF002248638452897 Merwe RVD Wan EA. The square-root unscented Kalman filter for state and parameter-estimation. in: Proceedings of (ICASSP '01) 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 6. Salt Lake City (2001). p. 3461–4. Julier SJ Uhlmann JK. Unscented filtering and nonlinear estimation. Proc IEEE (2004) 92:401–22. 10.1109/JPROC.2003.823141 Kalman RE. A new approach to linear filtering and prediction problems. Trans ASME J Basic Eng. (1960) 82:35–45. 10.1115/1.3662552 Silva FHLD Hoeks A Smits H Zetterberg LH. Model of brain rhythmic activity. Kybernetic (1974) 15:27–37. 10.1007/BF00270757 Faugeras O Touboul J Cessac B. A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs. Front Comput Neurosci. (2009) 3:1. 10.3389/neuro.10.001.200919255631 Grimbert F Faugeras O. Analysis of Jansen's Model of a Single Cortical Column. Research Report INRIA (2006). Pons AJ Cantero JL Atienza M Garcia-Ojalvo J. Relating structural and functional anomalous connectivity in the aging brain via neural mass modeling. Neuroimage (2010) 52:848–61. 10.1016/j.neuroimage.2009.12.10520056154 Silva FLD. 5. In: EEG: Origin and Measurement. World Scientific Publishing Co. (2011). p. 63–82. Available online at: http://www.worldscientific.com/doi/abs/10.1142/9789814525350_0005 Ary JP Klein SA Fender DH. Location of sources of evoked scalp potentials: corrections for skull and scalp thicknesses. IEEE Trans Biomed Eng (1981) 28:447–52. 10.1109/TBME.1981.3248177287042 Jurcak V Tsuzuki D Dan I. 10/20, 10/10, and 10/5 systems revisited: their validity as relative head-surface-based positioning systems. Neuroimage (2007) 34:1600–11. 10.1016/j.neuroimage.2006.09.02417207640 Berg P Scherg M. A fast method for forward computation of multiple-shell spherical head models. Electroencephalogr Clin Neurophysiol. (1994) 90:58–64. 10.1016/0013-4694(94)90113-97509274 Merwe RVD Wan EA. The unscented kalman filter for nonlinear estimation. In: Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000, AS-SPCC. Lake Louise, CA: IEEE (2000). p. 153–8. Solonen A Hakkarainen J Ilin A Abbas M Bibov A. Estimating model error covariance matrix parameters in extended Kalman filtering. Nonlinear Process Geophys. (2014) 21:919–27. 10.5194/npg-21-919-2014 Cantero JL Atienza M Gomez-Herrero G Cruz-Vadell A Gil-Neciga E Rodriguez-Romero R . Functional integrity of thalamocortical circuits differentiates normal aging from mild cognitive impairment. Hum Brain Mapp. (2009) 30:3944–57. 10.1002/hbm.2081919449329 Richards JE. Recovering dipole sources from scalp-recorded event-related-potentials using component analysis: principal component analysis and independent component analysis. Int J Psychophysiol. (2004) 54:201–20. 10.1016/j.ijpsycho.2004.03.00915331212 Crispin W Gardiner. Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences. Berlin: Springer (2004). Garcia-Ojalvo J Sancho J. Noise in Spatially Extended Systems, 1999th Edn. New York, NY: Institute for Nonlinear Science, Springer (1999). Easycap. Easycap Equidistant Layouts. (2018). Available online at: https://www.easycap.de/wordpress/wp-content/uploads/2018/02/Easycap-Equidistant-Layouts.pdf Toral R Colet P. Stochastic Numerical Methods: An Introduction for Students and Scientists. John Wiley & Sons (2014). Liu X Gao Q. Parameter estimation and control for a neural mass model based on the unscented Kalman filter. Phys Rev E. (2013) 88:042905. 10.1103/PhysRevE.88.04290524229247 Schiff SJ. Neural control engineering: the emerging intersection between control theory and neuroscience. In: Computational Neuroscience. MIT Press (2012). Available online at: https://books.google.es/books?id=P9UvTQtnqKwC Moye MJ Diekman CO. Data assimilation methods for neuronal state and parameter estimation. J Math Neurosci. (2018) 8:11. 10.1186/s13408-018-0066-830094571 Appendix: the unscented Kalman filter (UKF) algorithm

UKF is a predictor-corrector algorithm that estimates the state and parameters at a given time step k in two phases. The first one predicts the state based solely on the dynamical information of the system, i.e., the model. The second incorporates a measurement with which to correct the first estimation. Table A1 presents the symbols used in this paper for the variables of the Kalman filter.

Variables of the Unscented Kalman Filter.

Parameter Description
x^ State estimate
x^- a priori state estimate
Σ Sigma points
X* Transformed sigma points
X Redrawn sigma points [52]
Υ Sigma points projected to measurement space
y− Estimated measurement
z Measurement
f Model of the system
H Observation model
n State size
λ Scaling parameter
α Primary scaling factor
β Secondary scaling factor
κ Tertiary scaling factor
Wm Weight vector for the mean
Wcov Weight vector for the covariance
P State covariance estimate
P− a priori state covariance estimate
Pyy Predicted measurement covariance
Pxy State-measurement cross-covariance
Q State error covariance*
R Measurement error covariance*
K Kalman gain

The first step of the algorithm involves computing the expectation of the state and of the state covariance at time instant k + 1, known as the a priori estimation. For this we use a numerical implementation (using Heun's solver) of Jansen and Rit's model of a cortical column [26, 51], as described in the section 2.2.

The nature of the nonlinearities of this model prevents us from using a simple linearisation approach to propagating the statistics of the state variables across the transformation, as would be the case if we used the extended Kalman filter, for example. Therefore, we incorporate the unscented transform (UT) in our formulation of the Kalman filter, which, instead of attempting to propagate a distribution through the nonlinearity, first propagates a series of deterministically chosen points through the nonlinearity and then recovers the statistical information of the distribution from these.

Therefore, the a priori estimation of the state, x^k−, is obtained as follows, beginning with the calculation and projection of the 2n + 1 (where n is the state size) sigma points,

Σk−1,0=x^k−1 Σk−1,i=x^k−1+((n+λ)Pk−1)i,  i=1,...,n Σk−1,i=x^k−1−((n+λ)Pk−1)i−n,  i=n+1,...,2n

where Pk−1 is the estimated state covariance matrix for the previous time step. The square root of this matrix is well-defined, and can be calculated efficiently via a Cholesky decomposition [52]. This continues with the condensation of the projected sigma points into the a priori state estimate:

Xk|k−1*=f(Σk−1) x^k−=∑i=02LWim  Xi,k|k−1* Pk−=∑i=02LWicov[Xi,k|k−1*−x^k−][Xi,k|k−1*−x^k−]T+Q

where Q is the state error covariance and Wm and Wcov are the weight vectors, defined as

W0m=λn+λ W0cov=λn+λ+1−α2+β Wim=Wicov=12(n+λ),i=1,...,2n

In Equations 1 and 5, α, β and κ are scaling factors, and λ, which is crucial to guarantee a positive semi-definite covariance matrix P, is calculated as λ = α2(n + κ) − n. The primary scaling factor α determines the spread of the sigma points around the mean and is set at 0.001, it being usually set between 0.001 and 1 [63] and chosen according to the quality of the resulting estimation. The secondary scaling factor β contains prior information about the distribution of x; for Gaussian distributions, its optimal value is 2. Finally, κ, the tertiary scaling parameter, is set to 0, as is a usual practice [63].

We now use a measurement to correct the state estimation, which implies the mapping of the a priori estimate onto the measurement space for comparison. In our case, this transformation is a linear matrix H that relates the state of the cortical columns to an EEG reading (see section 2.2 for details). The sigma points Σk|k−1 are projected into the measurement space [52]

Υk|k−1=H[Σk|k−1],

from which the estimation of the measurement, y^k−, is calculated:

y^k−=∑i=02LWim  Υi,k|k−1

The second step of the algorithm corrects the a priori estimation of state and covariance by using the information available from the most recent measurement (in our case, an EEG reading). The impact of the measurement is determined by the Kalman gain Kk, which essentially expresses the level of confidence on the accuracy of the model and the level of noise in the data.

Pykyk=∑i=02LWicov  [Υi,k|k−1−y^k−][Υi,k|k−1−y^k−]T+R Pxkyk=∑i=02LWicov  [Xi,k|k−1−x^k−][Υi,k|k−1−y^k−]T Kk=Pxkyk  Pykyk−1 x^k=x^k−+Kk(zk−y^k−) Pk=Pk−−Kk  Pykyk  KkT

where Pykyk is the predicted measurement covariance, Pxkyk is the state-measurement cross-covariance, R is the measurement error covariance, and zk is the measurement for the current time step.

Funding. This work was partially supported by the Spanish Ministry of Economy and Competitiveness and FEDER (project FIS2015-66503) and by the Catalan Government (AGAUR grant FI-DGR 2014-2017). JG-O also acknowledges support from the Catalan Government (project 2014SGR0947), the ICREA Academia programme, and from the María de Maeztu programme for Units of Excellence in R&D (Spanish Ministry of Economy and Competitiveness, MDM-2014-0370).