Michelle Chong
email: chongms@unimelb.edu.au

Romain Postoyan
email: romain.postoyan@cran.uhp-nancy.fr

Dragan Nešić

Levin Kuhlmann
email: levink@unimelb.edu.au

Andrea Varsavsky

Estimating the unmeasured membrane potential of neuronal populations from the EEG using a class of deterministic nonlinear observers

Keywords: numbers: 87.19.le, 87.19.lj, 87.19.ll, 87.19.lr

We present a model-based estimation method to reconstruct the unmeasured membrane potential of neuronal populations from a single channel electroencephalographic (EEG) measurement. We consider a class of neural mass models that share a general structure, specifically the models by Stam et. al. in [1], Jansen et. al. in [2] and Wendling et. al. in [3]. Under idealised assumptions, we prove the global exponential convergence of our observer. Then, under more realistic assumptions, we investigate the robustness of our observer against model uncertainties and disturbances. Analytic proofs are provided for all results and our analyses are further illustrated via simulations.

Introduction

Brain studies involving the EEG typically employ signal processing strategies, which can be categorised into measure-based or model-based methods. The measure-based method involves extracting key measures from the EEG signal to classify brain phenomena. This approach is prevalent for the purpose of seizure detection (e.g. [START_REF] Osorio | Real-Time Automated Detection and Quantitative Analysis of Seizures and Short-Term Prediction of Clinical Onset[END_REF][START_REF] Gotman | Automatic detection of seizures and spikes[END_REF]) and prediction (e.g. [START_REF] Mirowski | Classification of patterns of EEG synchronization for seizure prediction[END_REF][START_REF] Kuhlmann | Patient-specific bivariate-synchrony-based seizure prediction for short prediction horizons[END_REF]) in epilepsy. These methods are continuously under development due to the strict conditions in which they can be considered successful [START_REF] Mormann | Seizure prediction: the long and winding road[END_REF].

On the other hand, model-based estimation methods can provide insight into the complex nature of the brain via mathematical models (e.g. [START_REF] Jansen | A neurophysiologically-based mathematical model of flash visual evoked potentials[END_REF][START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF][START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Babajani-Feremi | Multi-area neural mass modelling of EEG and MEG signals[END_REF] to name a few). These models describe brain activity at different scales, from a single neuron to populations of neurons and these activities are captured by states of a dynamical system. These states can be physiologically relevant and most are unmeasurable by conventional means. The model-based approach of estimating the unmeasured brain activity is called state estimation. In systems theory literature, the algorithm to estimate states from a measurement is commonly known as 'observer', 'filter' or 'estimator' [START_REF] Anderson | Optimal filtering[END_REF][START_REF] Kailath | Linear estimation[END_REF]. We will use the term 'observer' in this paper.

The objective of this study is to propose a model-based algorithm to estimate online the hidden states of neural mass models formulated as deterministic ordinary differential equations in continuous time. As such, this scheme similar in flavour to Dynamic Causal Modelling [START_REF] Garrido | Dynamic causal modelling of evoked potentials: A reproducibility study[END_REF][START_REF] Marreiros | Population dynamics: Variance and the sigmoid activation function[END_REF][START_REF] Marreiros | Population dynamics under the Laplace assumption[END_REF][START_REF] Moran | Dynamic causal models of steady-state responses[END_REF][START_REF] Daunizeau | Dynamic causal modelling: a critical review of the biophysical and statistical foundations[END_REF][START_REF] Marreiros | A dynamic causal model study of neuronal population dynamics[END_REF], which aims to recover the hidden states and parameters of an underlying dynamical system from observed time series. In this paper, we focus on estimating hidden states, noting that this is useful for model identification or parameter estimation. The ability to recover hidden states and parameters from electrophysiological data is of clear importance; for example, in the diagnosis, classification and anticipation of seizures in epilepsy.

Most of the work in state estimation for neuroscientific and neurological studies are developed under the stochastic framework, using the Kalman filter and its variants, or by employing the expectation maximisation algorithm [START_REF] David | Dynamic causal modeling of evoked responses in eeg and meg[END_REF][START_REF] Schiff | Kalman filter control of a model of spatiotemporal cortical dynamics[END_REF][START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF][START_REF] Ullah | Assimilating seizure dynamics[END_REF][START_REF] Lu | Polytopic linear parameter-varying model of epileptiform activity[END_REF][START_REF] Dr Freestone | A data-driven framework for neural field modeling[END_REF]. However, the usage of stochastic observers on a nonlinear model such as [START_REF] Jansen | A neurophysiologically-based mathematical model of flash visual evoked potentials[END_REF] and [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] requires that the observer be initialised close to the true initial condition, which is realistically unknown. Moreover, the convergence of the estimates to the true states is not guaranteed for every trajectory.

Deterministic state estimation circumvents these issues. Nevertheless, the main difficulty lies in that there exists no generic method of designing deterministic observers for nonlinear models. As existing techniques do not apply to the models we consider, we developed our own observer and provide analytical proofs for the convergence of state estimates to the true states.

Deterministic state estimation methods at the microscopic level (single neuron) include works by Tokuda et. al. in [START_REF] Tokuda | Parameter estimation for neuron models[END_REF], Totoki et. al. in [START_REF] Totoki | Firing pattern estimation of synaptically coupled hindmarsh-rose neurons by adaptive observer[END_REF], Tyukin et. al. in [START_REF] Tyukin | State and parameter estimation for canonic models of neural oscillators[END_REF] and Mao et. al. in [START_REF] Mao | Identification of biological neurons using adaptive observers[END_REF]. As our measurement of choice is the EEG, which is known to reflect the average behaviour of populations of neurons [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF], our work differs by taking a macroscopic view. We focus on designing deterministic observers for estimating the unmeasured mean membrane potential of neuronal populations. At present, it is only possible to measure the membrane potential of neuronal populations through invasive implantation of electrodes. The positioning and size of electrodes need to be accurate enough such that the measurements are not influenced by nearby populations. To this end, we consider physiologically relevant models that describe the interaction between neuronal populations, known as neural mass models (see [START_REF] Deco | The dynamic brain: From spiking neurons to neural masses and cortical fields[END_REF] for review), that are known to replicate patterns typically seen in the EEG (e.g. [START_REF] Jansen | A neurophysiologically-based mathematical model of flash visual evoked potentials[END_REF][START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF][START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF][START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF][START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] to name a few). We consider three neural mass models available in the literature that share a general mathematical structure, they describe: (i) alpha rhythms [START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF], (ii) alpha rhythms in the cerebral cortex [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] and (iii) epileptic activity in the hippocampus [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. We rewrite these neural mass models in state space form, specifically, in state coordinates that are amenable for observer design. We also identify the common model features that allow us to prove convergence of state estimates to the true states in a unified manner.

In the development of our observer, we first assume that there are no uncertainties and external disturbances that affect our model. We later relax these assumptions and show that the proposed class of observers is robust to uncertainty in parameters, additive disturbance, as well as noise and artefacts in the input and EEG measurement. The robustness of this class of observers to parameter uncertainty provides a starting point in the final aim of building an online observer for both the states and parameters of the model that will form the further work of this study. The contributions of this paper are:

• The design of a class of deterministic nonlinear observers for neural mass models that guarantee convergence of the estimated states to the true states for any initial condition and any known input. Analytical proofs using system theoretic methods are provided in the Appendix. • Robustness analysis of our proposed observers towards common issues faced in estimating brain activity. These include measurement noise, uncertainty in the parameters and the measured input. Simulation studies confirm the analyses performed.

The rest of the paper is organised as follows. We start by introducing the notations and formally state the problem in Section 2. Section 3 contains a more detailed exposition on the three neural mass models of interest. Then in Section 4, we show how to convert these neural mass models into a state space form that is amenable for observer design and error convergence analysis. Section 5 contains the main simplifying assumptions. Section 6 contains the main result of the paper that states the exponential convergence properties of the state estimates to the true states. We then relax some of the earlier strong assumptions and show that the observer is robust under a more realistic set of assumptions. We present simulation results in Section 7. In Section 8, We provide some discussions and conclude the paper.

Notation

• The matrix

      A 1 0 . . . 0 0 . . . . . . . . . . . . 0 0 . . . 0 A n       where A i for i ∈ {1, . . . , n} are m × m matrices, is denoted as diag(A 1 , . . . , A n ). • The set L ∞ denotes the set of functions f : R → R n , for some n ∈ Z, such that for any 0 ≤ t 1 ≤ t 2 < ∞ there exists r ≥ 0 so that f [t1,t2] := ess sup τ ∈[t1,t2] |f (τ )| ≤ r.
• A vector u ∈ R n , where each element is drawn from a Gaussian distribution with mean µ ∈ R and variance σ 2 ∈ R is denoted u ∼ N (M, σ 2 I), where M = (µ, . . . , µ) ∈ R n and I ∈ R n×n is an identity matrix.

Problem statement

We consider several classes of neural mass models described in [START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF], [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] and [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. These models have some common features and can all be represented by the same general functional block diagram that is shown in Figure 1. The input to the model u represents all the external influences from afferent neuronal populations or structures.

The output of the model y is taken to be the EEG measurement, which is proportional to the aggregated average postsynaptic potentials of the pyramidal neurons. The state x is a vector whose entries are the afferent population specific contributions to the mean membrane potential of each neuronal population and its derivative. The parameter θ is a vector whose entries are the connectivity strengths and the synaptic gain of each neuronal population. The state x, output y and parameter θ have subtly different meaning for each model we consider and these differences are discussed in more detail in the following sections. Our goal is to design a model-based algorithm (observer) which uses measurements from u and y to compute the state estimate x as illustrated in Figure 2. We show how to design the observer for a general class of systems in state space form that covers many neural mass models. The observer is designed to provide a state estimate x that converges to the real state x in exponential time, for any initial condition and any input. In the next section, we describe the neural mass models considered in more detail.

Neural mass models

We consider several classes of neural mass models that describe the interconnection between neuronal populations as shown in Figures 345. These block diagrams describe a functional relationship between various neuronal populations. They are developed by: (i) Stam et. al. in [START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF], (ii) Jansen and Rit in [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] and (iii) Wendling et. al. in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. These models are based on the models by Wilson and Cowan in [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF] as well as Lopes da Silva in [START_REF] Lopes | Models of neuronal populations: the basic mechanisms of rhythmicity[END_REF]. In the following sections, the classes of neural mass models are introduced with increasing level of complexity. Noting that the models by Stam We begin with the model by Stam et. al., which includes an excitatory and inhibitory population to replicate alpha rhythms as seen in Figure 3. The emergence of alpha rhythms in the EEG is related to the human subject being in a relaxed state with eyes closed. Our observer is used to estimate the membrane potential of each population and its derivatives from the EEG and input measurements. Estimating the unmeasured membrane potential of neuronal populations when alpha rhythms are seen on the EEG measurement could allow for the better understanding of the visual pathway while in an idle state.

Neural mass model by Jansen and Rit

The interconnections between the pyramidal neurons, excitatory and inhibitory populations are described in this model to investigate the generation of visual evoked potentials in the cerebral cortex, as illustrated in Figure 4. We design observers to estimate the unmeasured membrane potential of neuronal populations, which could give insight into the stimulation of a sensory pathway. 5. They are the pyramidal neurons, the excitatory population (included in the pyramidal neurons), the slow and fast inhibitory populations. The fast somatic projection of the inhibitory population is introduced in this model as it is hypothesised to play a role in the fast oscillatory pattern seen in the EEG during epileptic seizure onsets. Estimating the membrane potential contribution of one population to another could provide insight into the underlying dynamics of seizure activity.

!"#$%&'$() *+,#-./) 0.1&2&3-#") 0.3+#.+,#-./) Figure 6 describes the interaction between populations of neurons in greater detail, which consists of postsynaptic membrane potential (PSP) kernels h e , h i and h g , sigmoid functions S : R → R and connectivity constants C 1 to C 7 . The firing rate of the afferent population is converted into an excitatory, slow or fast inhibitory postsynaptic membrane potential via the following kernels, for t ≥ 0:

?) ?) @) !"#$%&'(%))'(*+"&' 0.1&2&3-#") 0.3+#.+,#-./) 4A$/3)/-%$67) 8#-9+76-.:) @) @) B C) B D) B E) B F) B G) B H) B I)
C1 C2 C3 C4 C5 ) (t h e ) (t h e ) (t h i ) (t h g ) (t h e ) (t h e ) (t h i ] . [ ! ! ! S 11 x !" #" #" #" !" $%&"'(")$)" !" !"
• The excitatory population:

h e (t) = θ A at exp(-at). (1)
• The slow inhibitory population:

h i (t) = θ B bt exp(-bt). (2)
• The fast inhibitory population:

h g (t) = θ G gt exp(-gt). (3)
Parameters θ A , θ B and θ G in (1)-(3) correspond to the synaptic gains of the excitatory, slow and fast inhibitory populations respectively. These parameters characterise the observed pattern in the EEG. For example, the values of θ A , θ B and θ G that distinguish between seizure and non-seizure activities have been identified in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. Internal variables x 11 , . . . , x 71 are introduced as shown in Figures 3456. They describe the membrane potential contribution from one population to another. For example as seen in Figure 6, the mean membrane potential of the pyramidal neurons is x 11 -x 21 -x 31 , which reflects the membrane potential contribution from the excitatory, slow and fast inhibitory populations respectively. The mean membrane potential of a population is converted into the average firing rate of all the neurons in that population using a sigmoid function S:

S(z) = α 2 1 + exp -r 2 (z -V 2 ) , for z ∈ R, (4)
where α 2 is the maximum firing rate of the population, r 2 is the slope of the sigmoid and V 2 is the threshold of the population's mean membrane potential.

The neuronal populations are connected with connectivity strengths C 1 to C 7 , which represents the average number of synaptic contacts between the neuronal populations concerned.

Neural mass models in state space form

Our observer design is most conveniently carried out using the state space form of the neural mass models. However, some of the neural mass models of interest presented in Section 3 were in block diagram form (e.g. Figure 6) in the given references [START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF][START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. In [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] and [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF], state space forms were provided but not in the convenient state coordinates where the techniques we use for proving convergence of state estimates can be applied. Therefore, we illustrate how this can be done on the model by Wendling et. al., whose detailed block diagram can be found in Figure 6. The other models considered are special cases of the model by Wendling et. al. and hence, can easily be obtained from the derivation below.

We will show that all neural mass models from Section 3 can be written in the following state space form:

ẋ = Ax + G(θ)γ(Hx) + σ(u, Cx, θ), (5)
and the output of the model is:

y = Cx, (6)
where the state vector is

x ∈ R n , input is u ∈ R, output/EEG measurement is y ∈ R, parameter vector is θ ∈ R p , G : R p → R n×m , nonlinearity γ = (γ 1 , . . . , γ m )
with γ i : R → R for i ∈ {1, . . . , m} and nonlinearity σ = (σ 1 , . . . , σ n ) with σ i : R×R×R m → R for i ∈ {1, . . . , n}. The number of states n, number of parameters p and number of scalar nonlinear functions m differs for each model. These are defined in Sections 4.3-4.5.

Essential features

We highlight the following important essential features that our neural mass models possess. These are crucial for the applicability of the nonlinear observer proposed in Section 6 to the neural mass models we consider in this paper:

(i) The eigenvalues of the matrix A in (5) have strictly negative real parts.

(ii) The nonlinearity γ in ( 5) is globally Lipschitz [START_REF] Khalil | Nonlinear systems[END_REF]Theorem 3.2] i.e. there exists ρ > 0 such that |γ(u)-γ(v)| ≤ ρ|u-v| for u, v ∈ R. Moreover, it is also bounded: γ(z) ≤ α, for all z ∈ R. For instance, this nonlinearity γ can take the form of (4), which satisfies both properties. (iii) The nonlinearity σ in ( 5) is bounded, i.e. for i ∈ {1, . . . , n}, there exists M > 0 such that |σ i (z)| ≤ M , for all z ∈ R.

The subsystems x i from (10) for i ∈ {1, . . . , 7} that make up the full model ( 5) are interconnected in a unique manner as shown in Figure 6. This interconnection leads to a 'cascaded' structure of the error dynamics system that is crucial in the convergence proofs stated in Appendix A. This cascaded structure will be introduced in Appendix A.

Physiological interpretation

Physiologically, the first term in (5) implements the postsynaptic potential (PSP) kernels from (1), ( 2) and (3). This is effectively a convolution of the pre-synaptic firing rates arriving from other populations with the appropriate PSP response functions. These firing rates are modelled in the second and third term in (5) that incorporates the sigmoid firing rate function of the depolarisation of contributing populations. The second term, G(θ)γ(Hx), reflects the influence of all states except the membrane potential of pyramidal population Cx. While the third term, σ(u, Cx, θ), reflects the influence of the mean membrane potential of the pyramidal cells Cx and the exogenous input u.

In the following sections, we present the state space form for each model for ease of observer design in Section 6. Detailed derivations are first shown for the model by Wendling et. al., as it is the most complex model and then the subtle differences in derivations are described for the other models.

State space form for the model by Wendling et. al.

We write the Wendling et. al. model in state space form by introducing the state variables x i1 for i ∈ {1, . . . , 7} as the membrane potential contribution from one population to another and x i2 for i ∈ {1, . . . , 7} as its derivative. The states x i1 is introduced at the outputs of all the impulse responses h e , h i and h g blocks as shown in Figure 6. Recalling that the Laplace transform of the impulse responses h e , h i and h g (as described by ( 1), ( 2) and ( 3)) are second order transfer functions, by performing the inverse Laplace transform, each transfer function is represented by a second order ordinary differential equation (ODE). We show this transformation for h e from (1) as an example. Let the input to the h e block be ū and output be ȳ. We denote the Laplace transform of signal v as L(v). Hence, the Laplace transform of h e is:

L(h e (t)) = L(θ A at exp(-at)) = θ A a (s + a) 2 . (7)
Recalling that L(h e ) = L(ȳ) L(ū) , we obtain:

L(ȳ)s 2 + 2aL(ȳ)s + a 2 L(ȳ) = θ A aL(ū). (8)
By taking the inverse Laplace transform, we obtain a second order ODE as follows:

ÿ + 2a ẏ + a 2 ȳ = θ A aū. (9)
The x i2 states are defined as x i2 = ẋi1 for i ∈ {1, . . . , 7} to rewrite the second order ODEs as two first order ODEs for each impulse response block. We illustrate this for the h e (t) block in the fast inhibitory population, then the output of that block is ȳ = x 51 and the input is ū = C 3 S(x 11x 21x 31 ). Taking x 52 = ẋ51 , (9) can be written as two first order ODE as follows:

ẋ51 = x 52 ẋ52 = -2ax 52 -a 2 x 51 + θ A aC 3 S(x 11 -x 21 -x 31 ).
Hence, each impulse response h e , h i and h g will each introduce a first order ODE in the following general state space form by taking

x i = [x i1 , x i2 ] T for i ∈ {1, . . . , 7}: ẋi = A i x i + 0, ϑ i S(µ i ) + ϕ i T , (10)
where µ i is the input to the respective sigmoid functions,

A i = 0 1 -k i1 k i2 -(k i1 + k i2 ) , for i = {1, . . . , 7} with k 11 = k 41 = k 51 = k 61 = a, k 12 = k 42 = k 52 = k 62 = a, k 21 = k 71 = b, k 22 = k 72 = b, k 31 = g and k 32 = g. ϑ i and ϕ i are defined as such ϑ 1 = θ A aC 2 , ϑ 2 = θ B bC 4 , ϑ 3 = θ G gC 7 , ϑ 7 = θ B bC 6 , ϑ 4 = ϑ 5 = ϑ 6 = 0 and ϕ 1 = θ A au, ϕ 2 = ϕ 3 = ϕ 7 = 0, ϕ 4 = θ A aC 1 S(y), ϕ 5 = θ A aC 3 S(y), ϕ 6 = θ A aC 5 S(y).
Constants a, b and g are strictly positive. S is a sigmoid function described by [START_REF] Osorio | Real-Time Automated Detection and Quantitative Analysis of Seizures and Short-Term Prediction of Clinical Onset[END_REF]. All constants discussed in this section are summarised in Appendix C. The subsystems defined in [START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF] are put together to be written compactly in state space form ( 5)-( 6) for ease of observer design in Section 6.

We take the state vector in ( 5) and ( 6) to be x = [x 1 , . . . , x 7 ] T where x i for i = {1, . . . , 7} satisfy [START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF]. The states x 1 , x 2 and x 3 capture the membrane potential contribution and its derivative of the excitatory, slow and fast inhibitory populations to the pyramidal neurons respectively. The states x 4 , x 5 and x 6 capture the membrane potential contribution and its derivative of the pyramidal neurons to the excitatory, slow and fast inhibitory populations respectively. The output is y = x 11x 21x 31 . The specific matrices in ( 5) and ( 6) are denoted as:

• The parameter vector is θ = [θ A , θ B , θ G ] T , • The matrix A = diag(A 1 , . . . , A 7 ), • γ = [S, S, S] T from (4), (11)
• σ = (0, θ A au, 0, 0, 0, 0, 0, θ A aC 1 S(y), 0, θ A aC 3 S(y), 0, θ A aC 5 S(y), 0, 0), where S is described by ( 4),

• C = [ 1 0 -1 0 -1 0 0 0 0 0 0 0 0 0 ], (12)
• G =                         0 0 0 θ A aC 2 0 0 0 0 0 0 θ B bC 4 0 0 0 0 0 0 θ G gC 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 θ B bC 6 0                         (13)
and

• H =  
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -1 0   .

State space form for the model by Jansen and Rit in [2]

We write the model in state space form by taking the state vector in (5) to be x = [x 1 , x 2 , x 4 , x 5 ] T , where x i for i = {1, 2, 4, 5} satisfy [START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF]. States x 1 and x 2 are respectively, the membrane potential contribution and its derivative of the excitatory and inhibitory populations to the pyramidal neurons. States x 4 and x 5 capture the membrane potential contribution and its derivative of the pyramidal neurons to the excitatory and inhibitory populations respectively. The output is y = x 11x 21 . The specific matrices in ( 5) and ( 6) are denoted as:

• The parameter vector is θ = [θ A , θ B , C 1 , C 2 , C 3 , C 4 ] T , • A = diag(A 1 , A 2 , A 4 , A 5 ),
• γ = [S, S] T where S is defined in ( 4), ( 14)

• σ = [0, θ A au, 0, 0, 0, θ A aC 1 S(y), 0, θ A aC 3 S(y)] T , • C = [ 1 0 -1 0 0 0 0 0 ], (15)
• G =             0 0 θ A aC 2 0 0 0 0 θ B bC 4 0 0 0 0 0 0 0 0             (16)
and

• H = 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 .

State space form for the model by Stam et. al. in [1]

The model is written in state space form by taking the state vector in (5) as x = [x 1 , x 2 , x 5 ] T where x i for i = {1, 2, 5} satisfy [START_REF] Robinson | Neurophysical modeling of brain dynamics[END_REF]. State x 1 represents the contribution and its derivative of the excitatory population's activity to itself. States x 2 and x 5 represent the contribution and its derivative of the inhibitory population to the excitatory population and vice versa respectively. The output is y = x 11x 21 . This model differs from the models by Wendling et. al. and Jansen and Rit in the sense that the firing rate of the population is converted to postsynaptic potential via different kernels from ( 1) and ( 2), for t ≥ 0:

• Excitatory population:

h e (t) = θ A [exp(-a 1 t) -exp(-a 2 t)]. (17)
• Inhibitory population:

h i (t) = θ B [exp(-b 1 t) -exp(-b 2 t)]. (18)
As performed in Section 4.3, by taking Laplace transformations of ( 17) and ( 18) and taking the inverse Laplace transform, the kernels can be written as second order ODEs. They can then be rewritten as two first order ODEs by introducing extra state variables, x i2 for i ∈ {1, 2, 5}, in a similar fashion as in Section 4.3. Also, the sigmoid function that converts the postsynaptic potential to the firing rate of the population differs from (4) for the models by Wendling et. al. and Jansen and Rit, as follows:

S 1 (z) =    α 1 exp r 1 (z -V 1 ) z ≤ V 1 , α 1 2 -exp -r 1 (z -V 1 ) z > V 1 , (19)
where α 1 is the maximum firing rates of the population, r 1 is the slope of the sigmoid and V 1 is the threshold of the population's mean membrane potential. Nevertheless, rewriting ( 17)-( 19) into state space form does not differ from the derivation presented in Section 4. The specific matrices in [START_REF] Gotman | Automatic detection of seizures and spikes[END_REF] are denoted as: [START_REF] Daunizeau | Dynamic causal modelling: a critical review of the biophysical and statistical foundations[END_REF], which still satisfies item (ii) in Section 4.1.

• The parameter vector is θ = [C 3 , C 4 ] T , • A = diag(A 1 , A 2 , A 5 ), • γ = S 1 from
(20)

• σ = [0, θ A (a 2 -a 1 )u, 0, 0, 0, θ A (a 2 -a 1 )C 3 S 1 (y)] T , • C = [ 1 0 -1 0 0 0 ], (21)
• G = [0, 0, 0, θ B (b 2 -b 1 )C 4 , 0, 0] T (22) and • H = [ 0 0 0 0 1 0 ].

Assumptions

Our aim is to design a class of nonlinear observers for the neural mass models considered in Section 3, such that the state estimates converge to the true states in exponential time, for any initial condition and any input. We first present the idealised assumptions in Section 5.1 that are used in Section 6.1. We then state the relaxed assumptions in Section 5.2 and we show the robustness of the observers in Section 6.2 under these relaxed assumptions.

Assumptions under the ideal scenario

The convergence of the state estimates provided by our observers are proven in Section 6.1 under the idealised assumptions as illustrated in Figure 2 and stated as follows:

Assumption 1. The synaptic gain of each neuronal population θ A , θ B , θ G and the connectivity strengths C 1 , C 2 , C 3 , C 4 are parameters that are constant and known.

These parameters are typically slowly-varying during a particular brain activity, such that they can be considered constant over the time period observed [START_REF] Khalil | Nonlinear systems[END_REF]Section 9.6]. When the brain transitions from one activity to another however, the parameters C 3 , C 4 for the model by Stam et. al. in [START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF], parameters θ A , θ B , C 1 , C 2 , C 3 and C 4 in for the model by Jansen and Rit in [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] and parameters θ A , θ B and θ G for the model by Wendling et. al. in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] will change. If the goal is solely to estimate the parameters, system identification methods such as least-squares estimation may be used [START_REF] Ljung | System identification: theory for the user[END_REF]. However, it remains a hard problem due to the nonlinearity of these models. For this study, we only consider the case where the brain is in one particular brain activity and that the parameters are known, that is for the case where the parameters have been identified a priori. We analyse the robustness of the observer to parameter uncertainty in Section 6.2.

Assumption 2. The input u is measured.

The input from afferent populations is hard to quantify in practice, hence this assumption is not justified in general. In Section 6.2, we relax this assumption by allowing the input to be uncertain. We show that good estimates can still be obtained provided that the L ∞ norm of the difference between the true and assumed input is small. This is further illustrated with simulations in Section 7.2.2. It is important to note that it is not an easy task to guarantee exponential convergence of the error estimates without assuming that the input is known exactly. Please refer to [START_REF] Hammouri | Unknown input observer for state affine systems: A necessary and sufficient condition[END_REF] and references therein. Assumption 3. The measured EEG y is noise-free.

The measured EEG considered here is recorded using intracranial electrodes, which we first assume to be clean data such that it can be considered noise-free. We will show later in Section 6.2 that the estimation error remains bounded despite noisy measurements. Assuming noise-free measurements makes it easier to design a deterministic observer and it is the first step to explain our methodology; we later relax our assumptions and show that the designed observer is robust to measurement noise. All the neuronal populations except for the pyramidal neurons or the excitatory population receiving the input u included in the models in Section 3 are not affected by adjacent populations. In Section 6.2, we will introduce the external influences as an additive disturbance and show that the estimation error remains bounded with bounded disturbance.

Relaxed assumptions

The assumptions set out in Section 5.1 are strong, but provide us with a good first step in proving the estimation error convergence of the observers. We now relax these assumptions to characterise their effects in Section 6.2 on the convergence property of the class of observers we propose in Section 6. The relaxed conditions are captured in Figure 7 and we restate our relaxed assumptions as follows:

Assumption 5 (Relaxation of Assumption 1). The parameters θ A , θ B , θ G , C 1 , C 2 , C 3 and C 4 of the model ( 5) are known with error.

We characterise these uncertainties in parameters by introducing bounded, timevarying signal θ (t) to the parameters of the model [START_REF] Gotman | Automatic detection of seizures and spikes[END_REF] as shown in [START_REF] Dr Freestone | A data-driven framework for neural field modeling[END_REF], so that the true parameter values in the model become θ + θ .

Assumption 6 (Relaxation of Assumption 2). The input is uncertain.

We introduce bounded, time-varying disturbances u (t) to the input available to the observer as shown in [START_REF] Tokuda | Parameter estimation for neuron models[END_REF]. Thus, the real input to (5) is modelled by u + u , whereas the observer is fed by an assumed input u. The introduction of u relaxes Assumption 2 by allowing the real input to be unknown to the observer.

Main results

We first propose a class of nonlinear observers in Theorem 1 under the ideal Assumptions 1-4. The ideal assumptions are later relaxed as stated in Assumptions 5-8. Theorem 2 uses the relaxed assumptions and show the robustness of the designed observers against disturbances and uncertainties.

Observer design under the ideal scenario

We propose the following class of observers:

ẋ = Ax + G(θ)γ(H x + K(C x -y)) + L(C x -y) + σ(u, y, θ), (23)
where x ∈ R n is the state estimate and observer matrices K ∈ R m , L ∈ R n are introduced via linear output injection terms K(C xy) and L(C xy).

The observer matrices L and K in [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] are the choice of the user. They are introduced to provide flexibility over the convergence rate of the estimates. We establish in Theorem 1 that for sufficiently small norms of matrices K and L, we obtain convergence of the estimates to the true states in exponential time, for any initial conditions. Theorem 1. Consider the model in general form [START_REF] Gotman | Automatic detection of seizures and spikes[END_REF] under Assumptions 1 to 4 and observer matrices K and L for ( 23) are chosen such that:

ρ|K||G| + |L| < ν 1 |C| , ( 24
)
where ρ is the Lipschitz constant of γ, matrices G and C are from the model ( 5) and ν > 0 is constructed in the proof. Then, the observer ( 23) is a global exponential observer for the model (5), i.e. for any u ∈ L ∞ and denoting the estimation error as e := xx,

|e(t)| ≤ k exp(-λt)|e(0)| ∀t ≥ 0, ∀e(0) ∈ R n , (25)
where constants k, λ > 0.

Proof. The proof is given in Appendix A.

Remark 1. Note that Theorem 1 applies to all neural mass models considered in Section 3, where matrices G and C are different for each model considered.

Consequently, the obtained observer matrices K and L as well as ν from ( 24) and k, λ from (25) will differ for each model.

The exponential convergence property of the observer is desirable in practice because it means that the estimates are guaranteed to converge to the true states in exponential time. Additionally, the validity of the convergence property for any initial conditions is highly desirable because the initial conditions of the neuronal populations are usually unknown.

Theorem 1 provides us with a condition (24) that gives a class of nonlinear observers parameterised by observer matrices K and L. This condition ( 24) is conservatively obtained in Appendix A. However, we see that the choice of L = 0 and K = 0 fulfils condition [START_REF] Ullah | Assimilating seizure dynamics[END_REF] and therefore, admits an open-loop observer for the model ( 5) [START_REF] Chong | A nonlinear estimator for the activity of neuronal populations in the hippocampus[END_REF]. The drawback of an open-loop observer lies in that the convergence speed of the state estimation error cannot be controlled by the user. Nevertheless, it remains a useful observer provided that the error converges sufficiently fast. In Section 7.1, we see in simulations that the open-loop observer for the model by Wendling et. al. [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] provides estimates that converge to the true states in a reasonable timeframe. Non-zero choices of K and L fulfilling condition [START_REF] Ullah | Assimilating seizure dynamics[END_REF] provide tuneability of the estimation error's convergence speed. However, condition [START_REF] Ullah | Assimilating seizure dynamics[END_REF] does not provide a priori information about the convergence rate of the observer. In Section 7.1, we show in simulations that some choices of non-zero matrices K and L that satisfy condition [START_REF] Ullah | Assimilating seizure dynamics[END_REF] can lead to faster convergence rate for the state estimation error.

Remark 2. The observer structure [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] shares the same mathematical structure as other nonlinear observers in the literature, that is, the high gain [START_REF] Raghavan | Observer design for a class of nonlinear systems[END_REF] (with K = 0) and circle criterion observers [START_REF] Arcak | Nonlinear observers: a circle criterion design and robustness analysis[END_REF]. These observers employ special techniques in obtaining observer matrices L and K that are not satisfied by the model we consider [START_REF] Gotman | Automatic detection of seizures and spikes[END_REF]. Therefore, we prove the existence of K and L matrices using a different analysis from that in [START_REF] Raghavan | Observer design for a class of nonlinear systems[END_REF][START_REF] Arcak | Nonlinear observers: a circle criterion design and robustness analysis[END_REF].

Robustness Analysis

We show in Theorem 2 that under Assumptions 5-8, i.e. uncertainties in modelling, parameters, measurement and input, the estimates provided are close to the true states, where the 'closeness' is determined by the L ∞ norm of the uncertainties y , θ , u and sys . The introduction of the uncertainties and disturbances to the ideal setup is illustrated in Figure 7.

Therefore, we obtain the following perturbed systems ( 26) and ( 27) from ( 5) and [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] as illustrated in Figure 7: From (5):

ẋ = Ax + G(θ + θ )γ(Hx) + σ(u, Cx, θ + θ ) + sys , y = Cx. (26)
From ( 23):

ẋ = Ax + G(θ)γ H x + K(C x -(y + y )) + L(C x -(y + y )) + σ(u + u , y + y , θ). ( 27
)
Theorem 2. Consider the perturbed estimation error system ( 26) and ( 27 This guarantees that for all e(0) ∈ R n , t ≥ 0 and for all u, y , θ , u and sys ∈ L ∞ , the error system satisfies the following:

|e(t)| ≤ k exp(-λt)|e(0)| + γ y ( y [0,t] ) + γ θ ( θ [0,t] ) + γ u ( u [0,t] ) + γ sys ( sys [0,t] ),
where k, λ > 0 and γ y , γ θ , γ u , γ sys : R → R are continuous positive increasing functions that are zero at the origin.

Proof. The proof for Theorem 2 is provided in Appendix B.

We see that with no uncertainties and disturbance, i.e. u = 0, θ = 0 and y = 0, we recover the results of Theorem 1. Loosely speaking, Theorem 2 states that small uncertainties and disturbance implies small estimation error. Remark 3. By setting u = -u, we have that the input to the observer ( 23) is 0, that is we consider the proposed observer [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] as an unknown-input observer. The estimation error of the proposed observer [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] converges with some error that depends upon the L ∞ norm of u . This result is advantageous because the mean firing rate of the afferent neuronal populations to a particular brain region is hard to measure in reality. We will further demonstrate this observation in Section 7.2.2 via simulations.

Simulation Results

We illustrate the performance of our observers using simulations. First, we perform simulations to show the convergence of the estimates under the ideal conditions stated in Section 5 and test the influence of the observation gains on the speed of convergence of the estimation error. Next, we show the robustness of our observers against uncertainties in modelling, parameters, input and measurement in Section 7.2.1. We also show that our observers can provide estimates that are close to to the true states when the input to the observer is unknown in Section 7.2.2. The simulation results provided are for the model of the hippocampus by Wendling et. al. in [3]. Similar results can be obtained for other models introduced in Section 3.

We select synaptic gains θ A = 5, θ B = 25 and θ G = 10 that correspond to seizure activity as identified in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. The initial condition of the model ( 5) is set to x(0) = [6, 0.5, 6, 0.5, 6, 0.5, 6, 0.5, 6, 0.5, 6, 0.5, 6, 0.5] T . We initialise the observer [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF] to x(0) = 0. This choice is arbitrary as we have shown that the convergence of the estimation error is valid for all initial conditions in Theorems 1 and 2.

In this section, we denote the observer matrices as

K = k × I k and L = l × I l where I k = [1, 1, 1] T and I l = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] T .
We call k and l the observer gains. For each of the scenarios in the following sections, the performance of two observers are evaluated:

(i) The open-loop observer l = 0, k = 0.

(ii) The observer with observer gains k = 0.1 and l = -0.2.

All other constants used in simulation are as specified in Appendix C.

Simulations under Assumptions 1 to 4: ideal scenario

We consider the ideal case where the parameters are constant and known, input u and measurement y are known and unperturbed as well as no modelling errors. Both model ( 5) and observer ( 23) are supplied with the same Gaussian noise input with mean 90mV and variance 30mV used in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF].

As stated in Theorem 1 and illustrated in Figures 8 and9, the state estimation error e := xx converges to 0 asymptotically in exponential time. We observe in simulations that for both observers, the convergence rate is in general faster than the duration of a specific brain activity (see [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]): by t = 0.3s, all state estimation errors have converged to 0. As seen in Figure 8, the plot of the error norm shows a large overshoot initially (approximately 1200). However, the plots on Figures 9 show that the estimates of the membrane potential contributions from one population to another xi1 are reasonably close to the true values x i1 . The large initial error is due to the estimation error in the other states, i.e. the error of the time derivative of the membrane potential contributions |x i2 -xi2 |, which are not the physical quantities of concern in practice.

As discussed in Section 6, the choice of non-zero observer gains provides flexibility over the convergence rate of the estimation error. Figure 8 shows that for a selection of observer gains k and l, the norm of the estimation error converges faster than with the case of k = l = 0 (open-loop observer). Next, we test the performance of the observers under uncertainties in the parameters θ A , θ B , and θ G , disturbances in the input and output of the observer, as well as additive disturbance.

We simulate the perturbed systems ( 26) and ( 27) with the perturbations: constant parameter uncertainty θ = (0.5, 2.5, 1), Gaussian input uncertainty u ∼ N (0, 0.3 2 ), Gaussian measurement noise y ∼ N (0, 0.1 2 ) and Gaussian model uncertainty sys ∼ N (0, 1 2 I). The performance of both observers are similar. Figures 10 and11 show that the estimation error converges to a neighbourhood of the origin. This observation agrees with Theorem 2.

Figure 10 shows that the estimates for the states of interest x i1 do converge reasonably close to the true states. It can be seen in Figure 11 that the relative error of the states is quite low (under 10%).

|x i -x i | max(x i )-min(x i ) for i ∈ {1, . . . , 7}
under the practical scenario (Assumptions 5 to 8) for the open-loop observer (grey solid line) and the observer with k = 0.1, l = -0.5 (black dashed line).

Simulations under Assumptions 1, 3, 4 and 6: uncertain-input observer.

The input to the brain is usually hard to quantify under realistic conditions. To gauge the performance of the observers under this condition, we perform another set of simulations where the input is set to u = 0 for the observers with no uncertainty in parameters and measurement noise. This case is formally stated in Remark 3.

For both observers, Figures 13 to 12 illustrate that the estimation error converges to a neighbourhood of the origin as expected. In Figure 13, the estimation error for subsystems (x 21 , x 22 ), (x 31 , x 32 ), (x 41 , x 42 ), (x 51 , x 52 ), (x 61 , x 62 ) and (x 71 , x 72 ) converges to a reasonably small neighbourhood of the origin. However, the (x 11 , x 12 ) subsystem exhibits a much larger steady state estimation error than the other subsystems due to the input u directly affecting it. Nevertheless, our simulations in Figure 14 show that the states can be reasonably reconstructed when the input is unknown.

In Figure 12, the open-loop observer outperforms the observer with k = 0.1 and l = -0.5, in the sense that the open-loop observer converges to a smaller neighbourhood around the origin e = 0 (with |e| ≤ 35 for the open-loop observer and |e| ≤ 85 for the observer with l = -0.5 and k = 0.1). A decision regarding the tradeoff between the convergence speed of the state estimation error and the accuracy of the estimates needs to be made when employing these observers. In this simulation scenario, the error converges to a neighbourhood in a reasonable timeframe (t = 0.2s) for both observers. Hence, the open-loop observer is the observer of choice due to the steady-state accuracy of its estimates.

Discussion and conclusion

Using a single channel EEG signal as the measurement of brain activity, we have successfully developed a class of nonlinear observers that can reconstruct the mean membrane potential of each neuronal population for several classes of neural mass models. We rewrote the neural mass models considered in suitable state coordinates and identified the essential model features that allowed us to prove that the proposed class of observers gives estimates that converge exponentially to the true states of the model, as stated in Theorem 1. Moreover, the global convergence property of our proposed observers ensures that any initial estimated state will converge to the true state of the model for any input. This result is desirable because we do not know the initial state of the brain. Limitations of our proposed class of observers lie in the requirement that parameters and input to the observers need to be known, both of which are not known in a clinical or experimental setting. However, we have shown in Theorem 2 that our observer is robust to uncertainties in parameters as well as disturbances in the model, input and EEG measurement, in the sense that bounded uncertainty will result in bounded estimation error.

Given that it is difficult to measure the actual mean membrane potential of each neuronal population, the estimates obtained serve to provide some understanding of the underlying dynamics during a specific brain activity. This is beneficial for the general EEG-based studies of sensory, motor and cognitive processing and can lead to the development of clinical treatment methods for brain disorders. One such example is in the treatment of epilepsy, where more positive outcomes of therapy may be achieved by targeting specific seizure-causing mechanisms.

So far, we have validated the efficiency of our observer analytically and further illustrated the results in simulation. The applicability of our observer to real EEG data will depend heavily on how good an abstraction the models considered are of the brain region concerned. Provided that the models are a good description of the brain region of interest, the estimates of the unmeasured brain activity provided by the observers from real EEG measurement are then a good depiction of real activity in the brain. Given that the model by Wendling et. al. [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] considered in this paper has been validated against EEG recordings obtained from depth electrodes implanted in the hippocampus of a temporal lobe epilepsy patient, we expect the estimates provided by our observer to be a good depiction of the underlying dynamics of a patient's epileptic brain. This is achieved under the assumptions stated in Section 5 and to a satisfactory degree under the more practical assumptions stated in Section 6.2.

It is important to note that neural mass models are 'lumped parameter' models that describe the properties of populations of neurons, as opposed to a single neuron. In this construct, neurons are assumed to be spatially close, so that their interconnections are dense enough such that the group of neurons can be considered as one population with an aggregated activity [START_REF] Soltesz | Computational neuroscience in epilepsy[END_REF]Chapter 23]. Based on these assumptions, the model captures the temporal evolution of neuronal population dynamics and neglects the spatial interactions. As such, the interpretation of the estimates for the model describing epileptic seizures originating from the hippocampus [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] for example, applies only to the genesis of seizures in a local sense and not to the propagation of seizures from a focus. Understanding the seizure propagating mechanisms will require distributed models, also known as neural field models such as those in [START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF][START_REF] Robinson | Dynamics of large-scale brain activity in normal arousal states and epileptic seizures[END_REF][START_REF] Liley | A spatially continuous mean field theory of electrocortical activity[END_REF] (see [START_REF] Deco | The dynamic brain: From spiking neurons to neural masses and cortical fields[END_REF] for a review). These distributed models are infinite-dimensional models that would require infinite-dimensional observers for the estimation of states. Turning an infinite-dimensional model to a lower order model that is finite-dimensional would allow for the development of suboptimal finite-dimensional estimation methods, which are far easier to develop and practical to implement. Ultimately, the choice of the model rests in the aim of the endeavour.

An alternative to using spatially distributed or continuous models is to interconnect neural mass models, like those considered here, into spatial networks. Such an approach has been employed in dynamic causal modelling for EEG by coupling Jansen and Rit neural mass models [START_REF] David | Dynamic causal modeling of evoked responses in eeg and meg[END_REF]. Our deterministic observers apply to a single neural mass model, but in future research could be extended to apply to interconnected neural mass models, and therefore be applicable to the popular dynamic causal modelling framework.

One of the motivations behind this work is to provide early detection of focal seizures in epilepsy. It has been identified in [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] that varying the synaptic gains of each neuronal population θ A , θ B and θ G leads to different EEG patterns. Our work is a starting point in tracking the transition of the brain from non-seizure to seizure activity. To meet this objective, the simultaneous estimation of both states and seizure-causing parameters will need to be performed. While estimating the parameters alone should achieve the aim of detecting seizure onsets, the additional estimation of states provides the added advantage of peeking into the hidden dynamics of each neuronal population's activity. The success in designing a class of nonlinear observers for a selection of neural mass models raises hope in tackling the challenging problem of detecting and predicting epileptic seizures. The approach presented in this paper is also applicable in general to neuroscientific studies, where estimating the unmeasured, or in some cases unmeasurable, activities of brain regions, from an EEG measurement is beneficial.

Appendix A. Proof of Theorem 1

We will do all proofs for the model by Wendling et. al. described in Section 3.3, which is the most general model that encompasses the models by Stam et. al. and Jansen et. al. described in Sections 3.1 and 3.2 respectively. The proof for all other models can be done by first identifying the cascade structure of the observation error system, then showing that the subsystems satisfy certain properties to conclude the global exponential stability of the whole observation error system. This is performed in a similar fashion in the proof for the Wendling et. al. model that follows.

From ( 5), [START_REF] Frogerais | Model and identification in epilepsy: from neuronal population dynamics to EEG signals[END_REF], the dynamics of the state estimation error e := xx is:

ė = Ae + G(γ(Hx) -γ(H x + KCe)) + LCe. (A.1)
The main idea is to consider the estimation error system (A.1) as the nominal error system (A.1) with L = 0 and K = 0 perturbed by the terms Gγ(H x) -Gγ(H x + K(C xy)) + LCe. As such, in Lemma 1, we build a Lyapunov function W for the nominal error system. Next, using W as a candidate Lyapunov function for (A.1), we obtain a bound for observer matrices K and L such that e = 0 is global exponential stable (GES) [START_REF] Khalil | Nonlinear systems[END_REF]Definition 4.5

].

Lemma 1. There exists a continuously differentiable W : R n → R such that the following holds, for all e ∈ R n : where k 1 , k 2 , k 3 and k 4 are strictly positive constants.

k 1 |e| 2 ≤ W (e) ≤ k 2 |e| 2 , ∂W (
Proof. The nominal estimation error system (A.1) with L = 0 and K = 0 is ė = Ae + G(γ(Hx)γ(H x)). (A.3) First note that system (5) has solutions that are defined for all time due to the global Lipschitz property of function S and the fact that input u ∈ L ∞ [START_REF] Khalil | Nonlinear systems[END_REF]Theorem 3.2]. In view of the same arguments, the solutions of system (23) (with y and u from system (5)) are also well-defined and exist for all time. Consequently, we have that (A.3) have solutions that are defined for all time.

We can decompose system (A.3) into 7 subsystems Σ e1 to Σ e7 :

Σ e1 : ė1 = A 1 1 + φ 1 (x 41 ) -φ 1 (x 41 ) Σ e2 : ė2 = A 2 e 2 + φ 2 (x 51 ) -φ 2 (x 51 ) Σ e3 : ė3 = A 3 e 3 + φ 3 (x 61 , x 71 ) -φ 3 (x 61 , x71 ) Σ e4 : ė4 = A 4 e 4 Σ e5 : ė5 = A 5 e 5 Σ e6 : ė6 = A 6 e 6 Σ e7 : ė7 = A 7 e 7 + φ 7 (x 51 ) -φ 7 (x 51 ), (A.4)
where e i := ( e i1 , e i2 ) = (

x i1 -xi1 , x i2 -xi2 ) ∈ R 2 for i = {1, . . . , 7}. φ 1 (x 41 ) = 0, θ A aC 2 S 2 (x 41 ) , φ 2 (x 51 ) = 0, θ B bC 4 S 2 (x 51 ) , φ 3 (x 61 , x 71 ) = 0, θ G gC 7 S 2 (x 61 - x 71 ) , φ 7 (x 51 ) = (0, θ B bC 6 S 2 (x 51
)). Matrices A 1 , . . . A 7 are as defined in Section 4 and have eigenvalues with strictly negative real parts. We note that subsystem Σ e1 is in cascade with Σ e4 , subsystems Σ e2 and Σ e7 are in cascade with subsystem Σ e5 , and subsystem Σ e3 is in cascade with Σ e6 and Σ e7 . Using this cascade structure, we show that the overall system (A.3) is GES by constructing the desired Lyapunov function W from the Lyapunov functions of each subsystems V i .

We consider Lyapunov functions V 1 , . . . , V 7 for each subsystem Σ e1 to Σ e7 of the form:

V i = e T i P i e i for i ∈ {1, . . . , 7},
where P T i = P i > 0 satisfies the Lyapunov equation P i A i + A T i P i = -I. This is always possible as the eigenvalues of A i have strictly negative real parts, in view of [START_REF] Khalil | Nonlinear systems[END_REF]Theorem 4.6]. We will show that each subsystem Σ ei is input-to-state stable (ISS) using V i [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF].

For subsystems Σ e4 , Σ e5 and Σ e6 , taking the derivative of V i = e T i P i e i along solutions of Σ ei for i ∈ {4, 5, 6}, we obtain:

Vi ≤ -|e i | 2 ≤ - 1 2 |e i | 2 .
Next, we show that V 1 is an ISS-Lyapunov function [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF] for subsystem Σ e1 w.r.t e 4 . Taking the derivative of V 1 = e T 1 P 1 e 1 along the solutions of Σ e1 , we obtain:

V1 = e T 1 (P 1 A 1 + A T 1 P 1 )e 1 + 2e T 1 P 1 (φ 1 (x 41 ) -φ 1 (x 41 )).
Since the function S 2 is globally Lipschitz with constant ρ 2 from (4), we have that:

|φ 1 (x 41 ) -φ 1 (x 41 )| ≤ θ A aC 2 ρ 2 |x 41 -x41 | ≤ ρ e1 |e 4 |, where ρ e1 = θ A aC 2 ρ 2 . Therefore, V1 ≤ -|e 1 | 2 + 2|e 1 ||P 1 |ρ e1 |e 4 |.
Recalling that ξχ ≤ 1 2 ξ 2 + 1 2 χ 2 , for any ξ, χ ∈ R and letting ξ = 2|P 1 |ρ e1 |e 4 | and χ = |e 1 |, we obtain the following:

V1 ≤ -|e 1 | 2 + 1 2 |e 1 | 2 + 2|P 1 | 2 ρ 2 e1 |e 4 | 2 ≤ - 1 2 |e 1 | 2 + 2|P 1 | 2 ρ 2 e1 |e 4 | 2 .
With similar arguments for the remaining systems, along the solutions of Σ ei , for i ∈ {2, 3, 7}:

Vi ≤ - 1 2 |e i | 2 + γ i |µ i | 2 ,
where µ 1 = e 4 , µ 2 = e 5 , µ 3 = (e 6 , e 7 ), µ 7 = e 5 and

γ i = 2|P i | 2 ρ 2 ei where ρ e2 = θ B bC 4 ρ 2 , ρ e3 = θ G gC 7 ρ 2 and ρ e7 = θ B bC 6 ρ 2 .
The composite Lyapunov function for the overall system (A.4) W is constructed using the Lyapunov function for each subsystem V i [START_REF] Sontag | Changing supply functions in input/state stable systems[END_REF]. We consider the candidate Lyapunov function that is positive definite and radially unbounded:

W = a 1 V 1 + a 2 V 2 + a 3 V 3 + V 4 + V 5 + V 6 + a 7 V 7 , (A.5)
for a 1 , a 2 , a 3 , a 7 > 0, which are determined below. The derivative of (A.5) along the solution of the overall error system (A.4) is

Ẇ ≤ -a 1 1 2 |e 1 | 2 + a 1 γ 1 |e 4 | 2 -a 2 1 2 |e 2 | 2 + a 2 γ 2 |e 5 | 2 -a 3 1 2 |e 3 | 2 + a 3 γ 3 (|e 6 | 2 + |e 7 | 2 ) - 1 2 |e 4 | 2 - 1 2 |e 5 | 2 - 1 2 |e 6 | 2 -a 7 1 2 |e 7 | 2 + a 7 γ 7 |e 5 | 2 . By taking 0 < a 1 < 1 2γ1 , 0 < a 2 < 1 γ2 1 2 -a 7 γ 7 , 0 < a 3 < 1 2γ3 and a 7 > 2a 3 γ 3 , we obtain Ẇ ≤ - 1 8γ 1 |e 1 | 2 - 1 2γ 2 1 2 -γ 7 |e 2 | 2 - 1 8γ 3 |e 3 | 2 - 1 2 |e 4 | 2 - 1 2 |e 5 | 2 - 1 2 |e 6 | 2 - 1 4 |e 7 | 2 . (A.6)
By letting P = diag(a 1 P 1 , a 2 P 2 , a 3 P 3 , P 4 , P 5 , P 6 , a 7 P 7 ) as well as denoting λ min ( P ) and λ max ( P ) as the maximum and minimum eigenvalues of P respectively, we have shown that (A.2) is fulfilled with k 1 = λ min ( P ) > 0, k 2 = λ max ( P ) > 0, k 3 = min{ We will show that the solutions of the error system (B.1) is input-to-state stable (ISS) [START_REF] Sontag | Input to state stability: Basic concepts and results[END_REF] with respect to the uncertainties y , u , θ and sys . For this purpose, we use function W as defined in (A.5).

The derivative of W along the solutions of (B.1) is: V 2 , α 2 , r 2 Parameters for the sigmoid function. α 2 is the maximum firing rate. r 2 is the slope of the sigmoid and V 2 is the threshold of the population's mean membrane potential. V 2 , α 2 , r 2 Parameters for the sigmoid function. α 2 is the maximum firing rate. r 2 is the slope of the sigmoid and V 2 is the threshold of the population's mean membrane potential. See [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] for values corresponding to different brain activity Table C3. Standard constants used in the model by Wendling et. al. in [3]

uFigure 1 .Figure 2 .

 12 Figure 1. Block diagram of the neural mass model.

Figure 3 .

 3 Figure 3. Neural mass model by Stam et. al. in[START_REF] Stam | Dynamics of the human alpha rhythm: evidence for non-linearity[END_REF]

Figure 5 .

 5 Figure 5. Neural mass model by Wendling et. al. in [3].

Figure 6 .

 6 Figure 6. Detailed block diagram of the Wendling et. al. model. Reproduced from Figure 4 in [3].

Assumption 4 .

 4 No external influences from other neuronal populations.

  Specifically, γ is[START_REF] David | A neural mass model for MEG/EEG: coupling and neuronal dynamics[END_REF] for the model byWendling et. al.,[START_REF] Kailath | Linear estimation[END_REF] for the model byJansen and Rit and (20) for the model by Stam et. al. Also, G is (13) and C is (12) for the model by Wendling et. al., G is (16) and C is (15) for the model by Jansen and Rit, and G is (22) and C is (21) for the model by Stam et. al..

Figure 9 .

 9 Figure 9. Membrane potential contribution x i1 (grey solid line) and estimated membrane potential contribution xi1 (red solid line: l = k = 0, black dashed line: k = 0.1, l = -0.5) under the ideal scenario (Assumptions 1 to 4).

7. 2 .

 2 Simulations for the practical scenario 7.2.1. Simulations under Assumptions 5 to 8: uncertainty in parameters, input and measurement as well as additive disturbance.

Figure 10 .Figure 11 .

 1011 Figure 10. True membrane potential contribution x i1 (grey solid line) and estimated membrane potential contribution xi1 (red solid line: l = k = 0, black dashed line: k = 0.1, l = -0.5) under the practical scenario (Assumptions 5 to 8).

Figure 12 .Figure 13 .

 1213 Figure 12. Norm of state estimation error |e| for the open-loop observer (grey solid line) and the observer with k = 0.1, l = -0.5 (black dashed line) under Assumptions 1, 3, 4 and 6.

Figure 14 .

 14 Figure 14. True membrane potential contribution x i1 (grey solid line) and estimated membrane potential contribution xi1 (red solid line: k = l = 0, black dashed line: k = 0.1, l = -0.5) under Assumptions 1, 3, 4 and 6.

  e) ∂e ≤ k 4 |e|, and along solutions to (A.1) with L = 0 and K = 0: Ẇ (e) ≤ -k 3 |e| 2 , (A.2)

1 a1 , 1 a2- 1 1, 1 b2- 1 C 3 , C 4

 1111134 LC)e + Ψ(x, x) + ∂W ∂e Ψ (x, x, y , u , θ , sys ).From Theorem 1, there exists k3 , k4 > 0 such that: Ẇ ≤ -k3 |e| 2 + k4 |e||Ψ |.Using the fact that the function S 2 is globally Lipschitz with Lipschitz constant ρ 2 and S 2 (z) ≤ α 2 for any z ∈ R as defined in (4), we obtain:Ẇ ≤ -k3 |e| 2 + k4 |e| σ y (| y |) + σ θ (| θ |) + σ u (| u |) + | sys | where σ y (| y |) = |(θ A aC 2 ρ, θ B bC 4 ρ, θ G gC 7 ρ, θ A aC 1 ρ, θ A aC 3 ρ, θ A aC 5 ρ, θ B bC 6 ρ)| + |L| | y |, σ u (| u |) = θ A a| u | and σ θ (| θ |) = | (aC 2 α + a u [0,t] )| θ |, bC 4 α| θ |, gC 7 α| θ |, aC 1 α| θ |, aC 3 α| θ |, aC 5 α| θ |, bC 6 α| θ | |. Therefore, if |e| > 2 k4 k3 σ y (| y |) + σ u (| u |) + σ θ (| θ |) + | sys | , .2),[35, (4.49) of Theorem 4.19] is fulfilled and [35, (4.48) of Theorem 4.19] is satisfied with α 1 = k 1 and α 2 = k 2 , where k 1 and k 2 are from Lemma 1 in Appendix A. Therefore, we can conclude that the error system (B.1) is ISS with respect to y , u , θ and sys with gains γ y (s) σ u (s) and γ sys (s) = k2 k1 2 k4 k3.s, for s ≥ 0 and γ y (0) = γ θ (0) = γ u (0) = γ sys (0) = 0. Here, we have shown it for the model by Wendling et. al.[START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF], but similar arguments apply to all other models considered in Section 3. Appendix C. Values and description of the constants Parameter Description Standard value Average time constant in the excitatory feedback loop a 1 = 55, a 2 = 605 s b1 Average time constant in the inhibitory feedback loop b 1 = 27.5, b 2 = 55 s -1 V 1 , α 1 , r 1 Parameters for the sigmoid function. α 1 is the maximum firing rate. r 1 is the slope of the sigmoid and V 1 is the threshold of the population's mean membrane potential. Average number of synaptic contacts in the inhibitory feedback loop C 3 = 32, C 4 = 3 θ A , θ B Synaptic gain of the excitatory and inhibitory populations respectively θ A = 1.65, θ B = 32Table C1. Standard constants used in the model by Stam et. al. in [1]. Average time constant in the slow inhibitory feedback loop b = 50 s -1

1 C 1 , C 2 C 3 = 1 g

 11231 1 , r 2 = 0.56 mV -Average number of synaptic contacts in the excitatory feedback loop With C = 135, C 1 = C and C 2 = 0.8 C C 3 , C 4 Average number of synaptic contacts in the slow inhibitory feedback loop C 4 = 0.25 C θ A and θ B Synaptic gain of the excitatory and inhibitory populations respectively θ A = 3.25 and θ B = 22 Table C2. Standard constants used in the model by Jansen and Rit in [2]. Average time constant in the slow inhibitory feedback loop b = 50 s -1 Average time constant in the fast inhibitory feedback loop g = 500 s -1

- 1 C 1 , C 2 7

 1127 Average number of synaptic contacts in the excitatory feedback loopWith C = 135, C 1 = C and C 2 = 0.8 C C 3 , C 4Average number of synaptic contacts in the slow inhibitory feedback loopC 3 = C 4 = 0.25 C C 5 , C 6Average number of synaptic contacts between the fast and slow inhibitory feedback loop C 5 = 0.3 C and C 6 = 0.1 C C Average number of synaptic contacts in the fast inhibitory feedback loop C 7 = 0.8 C θ A , θ B and θ G Synaptic gain of the excitatory, fast inhibitory and slow inhibitory populations respectively

  et. al. as well as Jansen and Rit can be obtained from the model by Wendling et. al. by omitting certain neuronal populations, a detailed block diagram is presented only for the Wendling et. al. model in Section 3.3. 3.1. Neural mass model by Stam et. al.

$%&'(%))'(*+"&'

	+,7.%8* u(t)	=* !"#? @*	=* ? A*	0)(&1$2&3* 4-.(',/*	=* ? C*	>* ? B*	9.%7.%8* y(t) :!!;<*
		!"#$%&%'()*		+,5$6$%'()*
		+,%-(,-.(',/*		+,%-(,-.(',/*
	Figure 4. Neural mass model by Jansen and Rit in [2]
	3.3. Neural mass model by Wendling et. al.
	Wendling et. al. built upon the Jansen and Rit model described in Section 3.2.
	Four neuronal populations (with one population being a subset of another) are
	included in this model as shown in Figure

$%&"'()'(*&

				31%+,/"&,22"&'.1/"	)*+,&-.,/"01%+'02"
							41 x
	!"#$%&+,)'(*&	u	(t	31 x !"#) ] . [ ! ! ! S C7"	y	(t	)
						21 x
				71 x 405-6-7'+*"-071+01%+'02" 82/'9".10.+-:;" <+'=1;:'0>"	C6	S	] [ ! . ! !	51 x
								405-6-7'+*"-071+01%+'02"
				S	] [ ! . ! !	61 x		8(,27"2'&,:;"<+'=1;:'0>"

  Assumption 7 (Relaxation of Assumption 3). The measured EEG is affected by measurement noise.Disturbance y is introduced to characterise measurement noise.The measurement with noise is denoted as y + y .Assumption 8 (Relaxation of Assumption 4). The neural mass model is inaccurate.Additive disturbance sys (t) characterises the possible influence from other populations to the neuronal populations included in the model.

			Parameter uncertainty, ) (t	Model uncertainty, ) (t sys
						Measurement disturbance, ) (t
	Model input,	u		Neural mass model	Model output,	y
	Input		x	,		(noise--free EEG measurement)
	uncertainty, ) (t u
		State estimate,	x	Observer
						EEG measurement
						with noise

y Figure 7. Model and observer setup under the relaxed assumptions.

  ) under Assumptions 5 to 8. The observer[START_REF] Tokuda | Parameter estimation for neuron models[END_REF] with observer matrices K and L are chosen such that[START_REF] Ullah | Assimilating seizure dynamics[END_REF] is satisfied. Specifically, γ is (11) for the model by Wendling et. al., (14) for the model by Jansen and Rit and (20) for the model by Stam et. al. Also, G is (13) and C is (12) for the model by Wendling et. al., G is (16) and C is (15) for the model by Jansen and Rit, and G is (22) and C is (21) for the model by Stam et. al..

  } > 0 and k 4 = 2| P | > 0. κ 2 Ce , 0, θ G gC 7 S 2 (x 61 -x 71 )-S 2 (x 61 -x 71 -κ 3 Ce) , 0, 0, 0, 0, 0, 0, 0, θ B bC 6 S 2 (x 51 )-S 2 (x 51κ 2 Ce) andΨ = 0, θ A aC 2 S 2 (x 41κ 1 Ce) -S 2 (x 41κ 1 Cey ) + aC 2 S 2 (x 41 ) θ + au θ + θ A a u , 0, θ B bC 4 S 2 (x 51 -κ 2 Ce)-S 2 (x 51 -κ 2 Cey ) +bC 4 S 2 (x 51 ) θ , 0, θ G gC 7 S 2 (x 61 -x71 -κ 3 Ce)-S 2 (x 61 -x71 -κ 3 Cey ) +gC 7 S 2 (x 61 -x 71 ) θ , 0, θ A aC 1 S 2 (y)-S 2 (y + y ) + aC 1 S 2 (y) θ , θ A aC 3 S 2 (y) -S 2 (y + y ) + aC 3 S 2 (y) θ , 0, θ A aC 5 S 2 (y) -S 2 (y + y ) + aC 5 S 2 (y) θ , 0, θ B bC 6 S 2 (x 51κ 2 Ce) -S 2 (x 51κ 2 Cey ) + bC 6 S 2 (x 51 ) θ + sys .

	1 8γ1 , 1 2γ2	1 2 -γ 7 , 1 8γ3 , 1 4

Acknowledgements

This work was carried out with support from the St. Vincent's Hospital (Melbourne) Research Endowment Fund.

Continuing the proof of Theorem 1:

The derivative of W along the solutions to (A.1) is:

From Lemma 1:

As

where k 3 and k 4 are constructed in Lemma 1. Therefore, if K and L satisfy the following condition:

Therefore, the origin of the estimation error system (A.1) is GES according to [35, Definition 4.5], i.e. for all t ≥ 0,

for k, λ > 0.

Appendix B. Proof of Theorem 2

The proof that follows is performed for the Wendling at. al. model, where the models by Stam et. al. as well as Jansen and Rit can be derived from. The proof for these models can be performed in a similar fashion. From ( 26) and ( 27), the perturbed error system is: