Analysis and Control of Epileptiform Spikes in a Class of Neural Mass Models

The problemof analyzing and controlling epileptiform spikes in a class of neuralmassmodels is concernedwith. Since themeasured signals are always contaminated by measurement noise, an algebraic estimation method is utilized to observe the state from the noisy measurement. The feedback control is constructed via the estimated state. The feasibility of using such a strategy to control epileptiform spikes in a regular network of coupled neural populations is demonstrated by simulations. In addition, the influence of the type of the controlled populations, the number of the controlled populations, and the control gain is investigated in details.


Introduction
Since the 1970s, the use of neural mass models has been the preferred approach to model magnetoencephalography (MEG) and electroencephalography (EEG) signals.These models consist of macrocolumns or even cortical areas, describing the average activity of the whole population with a small number of state variables, which summarize the behavior of millions of interacting neurons [1].The original idea is put forward by Lopes Da Silva et al. [2] to use these models for alpha rhythms generation.Starting from Lopes Da Silva's model, Jansen et al. design a computer model that generates spontaneous EEG and evoked potentials [3].Subsequently, they extend the model in order to produce EEG and visual evoked potentials from two coupled populations of neurons in the visual cortex [4].Based on Jansen's model, neural mass models with different forms are developed to account for more complex transient and oscillatory behavior [1,5,6], to investigate effective connectivity from EEG or fMRI data [7], and to study the relations between brain rhythm changes and the connectivity among the regions of interest (ROIs) [8,9], or how event-related dynamics depend on extrinsic connectivity [10].Recently, dynamical causal modeling (DCM) [11][12][13] has been developed for the analysis of event-related [14,15] and steady-state responses [16].
A substantial progress in the use of Jansen's model is to simulate the dynamics of real EEG signals measured in the hippocampus during the transition from interictal to fast ictal activity of epilepsy seizures [17].An epileptic neural mass exhibits features of excitability which is reflected in its use to produce interictal spiking activities in a single population in response to random input [17].Studies also show that spiking activity can be propagated between populations in networks of two and three coupled neural populations [17].An extended neural mass model containing two explicitly modeling time scales of inhibition is provided by Goodfellow et al. to produce slow spike wave discharges (SWD) as well as sinusoidal background oscillations, thus providing the method for a parameter-driven transition to epilepsy seizure [18].The effect of local functional heterogeneities which are considered crucial for the mechanisms of epilepsy [19] is investigated in a spatially extended neural mass formulation with local excitatory and inhibitory circuits [20].These investigations play important roles in understanding the mechanisms of epilepsy seizures.However, until now, clinical electrotherapy for seizures relies heavily on empirical tuning of parameters and protocols [21,22], the effects of which are usually not optimal.Recently, it has been reported that the closed-loop control via feedback can optimize the entire seizure control process including efficacy of treatment, minimization of side effects, improvement of response by providing intermittent or minimal stimulation, reduction of damage, and minimization of power consumption [23].In closed-loop control strategy, control effects are evidently state specific.Thus, the development of feedback algorithms requires precise metric of system state.The control tools developed and implemented in modern control engineering, such as modern aeronautics, robotics, and manufacturing plants, can be used to predict the state of neural systems and control epilepsy seizures [24].
Neural activity is always measured through observing just a single variable such as voltage.A combination of noise in neurons and amplifiers as well as uncertainties in recording equipment such as electrode resistance and capacitance results in uncertainty of the measurement.In this study, the feedback closed-loop control based on algebraic estimation is utilized to control epileptiform spikes in Jansen's neural mass models.The algorithm of algebraic estimation [25][26][27] plays a role of filtering, which observes the system state.Comparing with other filtering methods, the attraction of algebraic estimation is its simplicity and real time of implementation and may be model-free.For a network of coupled neural populations, it is possible to change dynamical behaviors by imposing the control action on different populations.Our goal here is to establish rules for proper control of suppressing epileptiform spikes from a regular network of coupled neural populations using the proposed control strategy.The key point is to find the more appropriate populations over which the control action should be imposed.Thus, the rules might be related to the type and the number of populations that are selected to receive the function of feedback control.These relations are investigated and reported here by using plentiful numerical simulations.

System Description.
The basic idea of Jansen's neural mass model is to make excitatory and inhibitory populations interact such that oscillations emerge.A cortical area is understood as the constituent of three different populations of neurons, that is, excitatory pyramidal cells, excitatory stellate cells, and inhibitory interneurons.The excitatory pyramidal cells receive inhibitory and excitatory feedback from intrinsic interneurons and excitatory input from extrinsic areas.Any extrinsic input is characterized by a pulse density () which is relevant to time .The evolution of the population dynamics relies on two operators.The first one transforms the average pulse density of action potentials into an average postsynaptic membrane potential with the impulse response function given by for excitatory case and for inhibitory case, where () is the Heaviside function,  and  are average excitatory and inhibitory synaptic gains which determine the maximal amplitude of the postsynaptic potentials, and  and  are connected to both the membrane average time constant and the average distributed time delays in the dendritic tree.If time constants 1/ and 1/ are fixed,  and  can be applied to adjust the sensitivity of excitatory and inhibitory synaptic, respectively.Equations ( 1) and ( 2) can be represented by second-order ordinary differential equations where (), () are the input and output signals.By introducing a new variable  1 () and letting ż () =  1 (), (3) can also be written as follows: The second operator transforms the average membrane potential of the population into an average pulse density of potentials fired by the neurons, which is described by the nonlinear sigmoid function where 2 0 is the maximum firing rate, V 0 is the postsynaptic potential corresponding to a firing rate  0 , and  is the steepness of the sigmoid.Interactions between pyramidal cells and interneurons are summarized by four connectivity constants  1 −  4 , which account for the intrinsic circuitry and the average number of synaptic contacts.The above description is related to single neural population and summarized in dashed block in Figure 1.
It should be noted that the single neural population model concentrates on the brain's dynamical behaviors of a separate area.The spatial characters of EEG signals cannot be simulated by this model.Multiple coupled neural populations model (Figure 2) is appropriate to model the spatial characters of EEG because it takes into account the interaction between two or more cortical areas through introducing parameters related to the ways of populations connected and the delays associated with these connections [1,17].In Figure 2, population  is described concretely by Figure 1.As is observed in Figures 1 and 2, a gain constant   is introduced to define the strength of coupling between population  and population , while a filter with an impulse response ℎ  () = () ⋅    −   is introduced to model the delay of connections from population , where   is the average time delay on efferent connection from population .Appropriate setting of parameters   allows the populations to be interconnected unidirectionally or bidirectionally.The function ℎ  () is able to be rewritten as the form of (4).Therefore, multiple coupled neural populations model is described  by a set of eight first-order ordinary differential equations per population where the subscript  ( = 1, 2, . . ., ) represents the population under consideration, population  receives afferent information from population  ( = 1, 2, . . ., ,  ̸ = ) and the neighborhood (i.e., the pulse density   ()), and the average depolarization of pyramidal cells  3 () −  5 () is the output of population  and is assumed to model the EEG signal.The standard value of the model parameters is given anatomically as [4] As described by Jansen and Rit [4], each population in (7) produces well-defined alpha activity under the standard parameters.However, a modification of the parameters can result in the generation of epileptiform spikes as shown in simulation section.

Algebraic Estimation of Noisy Signals.
The method of algebraic estimation of noisy signals is addressed.The key point of the method given here lies in solving a classical polynomial approximation of the signal.Consider an analytic signal (), which is not directly available but rather is observed through a noisy measurement.In the sequel,   () = () + () is set for the noisy observation, where   () denotes the noisy measurement, and () denotes some additive noise.
The aim is to estimate () and its derivative from   ().
The convergent Taylor series expansion of () at  = 0 is as follows: The signal () and its derivative are approximated by the truncated Taylor polynomial of order , that is, Using elementary operational calculus,   () is transformed into The key idea of algebraic estimation is to isolate the th coefficient  () (0),  = 0, 1, 2, . . .,  by applying an appropriate operator.As a matter of fact, In light of Leibniz' formula for the differentiation of products, ( 12) is turned into where The right side expression of  in (13) matches a respective convolutional integral.Hence, ( 13) is carried out to the backward transform to the time domain where  > 0 is the interval of integration, The calculation of   (0) is derived from an integral over the time interval [0, ] for a small .According to (  ( − ))/  | =0 = (−1)   () (),  () () is expressed as a convolution product Replacing   in (17) by   , an algebraic estimator [ () ()]  of  () () is derived from the noisy measurement For specific , , and V, different polynomials for estimating the signal are obtained, such as The value of the integral in ( 18), (19), and ( 20) may be approximated by a trapezoidal numerical integration.For each sample time  =   with  = 0, 1, 2, . .., it follows: where   is the sample time,  = /  denotes the number of summation steps, Π  = Π V (,   ).
A feedback controller based on algebraic estimation is carried out to control epileptiform spikes in the neural mass model, as shown in Figure 3, where () is the output of the neural mass model, () is the noise,   () is the noisy measurement, [()]  is the algebraic estimation of (),  is the feedback gain matrix, and () is the control law with the form The goal is to control epileptiform spikes with as less control energy as possible.In this context, the control energy is defined as where the superscript  denotes the transpose of a vector or matrix.

Analysis of Epileptiform Spikes in the Neural Mass Models.
Mechanisms of epileptic seizures have been studied widely in recent years.It has been reported that interictal or ictal spike discharges in focal or partial seizures originate in an area of the cortex that is excessively excitable [28], and the interactions between neural populations are found to be strengthened during epilepsy seizures [29].In the neural mass model, the ratio / controls the degree of excitatory within a given population and directly influences the type of signal dynamics produced by the model.The value of this ratio can be altered by adjusting .The connectivity constant   can be used to adjust the strength of coupling between population  and population .This constant may also be used to interconnect populations in different ways.Consider a regular network represented by (7) with three identical neural populations, as shown in Figure 4. Several simulations are provided to validate the ability of the model to produce epileptiform spikes and clarify the causality relations of the produced signals among underlying populations.In our simulations, ( 7) is solved by using a fourth-order Runge-Kutta differential solver, and   ()( = 1, 2, 3) is assumed to be a Gaussian white noise with mean value 101 and standard deviation 35.Population 1 is made hyperexcitable by keeping all parameters standard except for the average excitatory synaptic gain , which is set to 3.4 mV.Populations 2 and 3 exhibit normal activity by keeping all parameters standard.Figure 5 shows the simulated dynamics of the model with different values of coupling constants.When  12 =  23 =  31 = 0, population 1 produces sporadic spikes appearing randomly while populations 2 and 3 produce signals reflecting normal activities (Figure 5(a)).Thus, the increase of  may result in the production of sporadic spikes.When  12 =  23 = 100,  31 = 0, the propagation of sporadic spikes from population 1 to populations 2 and 3 is observed, as shown in Figure 5(b).It may be conjectured that sporadic spikes in populations 2 and 3 are due to that in population 1 and the increase of the coupling constants  12 and  23 .When  12 =  23 =  31 = 100, a simulated EEG series with spikes resembling real EEG signals during the propagation of temporal lobe seizure [17] is observed, as shown in Figure 5(c).The introduction of  31 is sufficient to generate sustained spike discharges.Based on these results, the conclusion is reached that the neural mass models can generate epileptiform EEG signals such as sporadic spikes and sustained discharge of spikes by increasing the value of  and introducing the interactions between populations.The coupling strength and the coupling direction control the causality relations of the produced signals among the underlying populations.This observed behavior is in accordance with results obtained in [17].

Control of Epileptiform Spikes in the Neural Mass models.
The feedback control based on algebraic estimation is implemented to control epileptiform spikes in the regular network of three coupled populations (Figure 4).For the underlying model given in Figure 4, we choose a specific control strategy with the feedback gain matrix  of the following forms: where  1 = ( 0 0 0  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ) ,  2 = ( 0 0 0 0 0 0 0 0 0 0 0  2 0 0 0 0 0 0 0 0 0 0 0 0 ) ,  3 = ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  3 0 0 0 0 1 ,  2 , and  3 are constants.The rules of determining   are to control epileptiform spikes with as less control energy as possible.Several numerical simulations are applied to determine the following characteristics regarding the proposed control strategy in underlying neural mass model: (1) the relationship between the type of populations to be controlled and the ability to suppress epileptiform spikes, (2) the relationship between the control energy and the number of controlled populations, and (3) the dependence of the control time and the control energy on the control gain.In all cases, () is chosen as a normal distribution of random noise, and the polynomial ( 19) is used to derive the estimation []  from   .The parameters  = 0.25 s, V = 0, and   = 0.0025 s are set.

Effect of the Type of Controlled Populations.
The effect of adding feedback to different types of populations on the ability of controlling epileptiform spikes is analyzed.The ability of controlling epileptiform spikes is evaluated based on the output waveforms of the model for two different types of populations to be controlled: (1) hyperexcitable population by keeping all parameters standard except for , which is set to 3.4 mV, and (2) normal population by keeping all parameters standard.The model that contains one (population 1), two (populations 1 and 2), and three hyperexcitable populations is taken into account for a thorough analysis.In these cases, the coupling constants  12 =  23 =  31 = 100 are fixed.
For the model containing one hyperexcitable population, two different control schemes are performed to realize the goal.The first one is to add control action to the hyperexcitable population 1 by setting  1 = 1.96,  2 =  3 = 0.The second one is to add control action to the normal population 2 by setting  2 = 6,  1 =  3 = 0. Simulation results are exhibited in Figure 6.It is illustrated from the output waveforms that the first scheme is able to control spikes in all populations, while the second one fails to control spikes in population 1.In our simulations, control action is also added to population 3 or both the normal populations; spikes in population 1 cannot be controlled.
For the model containing two hyperexcitable populations, control action is added to hyperexcitable population 1 by setting  1 = 5.5,  2 =  3 = 0, the normal population 3 by setting  1 =  2 = 0,  3 = 10, and both the hyperexcitable populations by setting  1 =  2 = 0.86,  3 = 0, respectively.By comparing the output waveforms (as shown in Figure 7) for the applications of the above three control schemes, it is illustrated that epileptiform spikes in all populations are controlled only by adding control action to both hyperexcitable populations.
For the model containing three hyperexcitable populations, control action is added to hyperexcitable populations 1 and 2 by setting  1 =  2 = 8,  3 = 0 and all hyperexcitable populations by setting  1 =  2 =  3 = 1.62.The output waveforms for the use of the above two control schemes are exhibited in Figure 8. Spikes in all populations are controlled by adding control action to all hyperexcitable populations.
As such, these simulation results indicate that for the analyzed model the hyperexcitable populations are the necessary choices for being controlled to realize the goal by using the specific control strategy with the feedback gain matrix given in (24) and (25).Nevertheless, it is important to use as less control energy as possible to realize suppression of epileptiform spikes.The control energy is associated with the number of populations to be controlled from the definition of it.This association is investigated in details as shown next.

Effect of the Number of Controlled Populations.
The effect of the number of controlled populations on the control energy is analyzed.For the entire analysis, the coupling constants  12 =  23 =  31 = 100 are set.
First, population 1 is assumed to be hyperexcitable by keeping all parameters standard except for , which is set to 3.4 mV, and populations 2 and 3 are assumed to exhibit normal activity by keeping all parameters standard.Based on the analysis of Section 3.2.1, it is known that under the control strategy with ( 24) and ( 25) adding feedback to all hyperexcitable populations is necessary to control epileptiform spikes.Thus, the following three schemes in which the control action works on different numbers of populations are performed for comparison: (1) adding control action to population 1 by setting  1 = 1.96,  2 =  3 = 0, (2) adding   control action to populations 1 and 2 by setting  1 = 0.8,  2 = 0.2,  3 = 0, and (3) adding control action to all populations by setting  1 = 0.7,  2 = 0.35, and  3 = 0.175.Table 1 presents the total control energy over 20 simulation realizations for the above four control schemes.It is noted that the total control energy decreases with the increase of the number of controlled populations.Then, populations 1 and 2 are assumed to be hyperexcitable by keeping all parameters standard except for , which is set to 3.4 mV, and population 3 is assumed to exhibit normal It is concluded from the above results that under the control strategy with the feedback gain matrix given in ( 24) and ( 25) the total control energy for controlling epileptiform spikes in the analyzed model reduces when the number of controlled populations increases.

Effect of the Control Gain.
The effect of increasing control gain values   ,  = 1, 2, 3 on the control time and the control energy is analyzed.It has been shown in Section 3.2.2 that adding control action to each population needs less control energy for the regular network of three coupled neural populations.Thus, the effect of control gain is analyzed under this control frame.Here, the control time is approximately evaluated based on the moment that spikes in the model disappear.
Population 1 is assumed to be hyperexcitable by keeping all parameters standard except for , which is set to 3.4 mV and populations 2 and 3 are assumed to exhibit normal activity by keeping all parameters standard.For the convenience of analysis, let  2 =  1 /2 and  3 =  1 /4. Figure 9 presents the output waveforms of population 2 for different values of control gain  1 .The output waveforms of populations 1 and 3 are omitted because they are similar to that of population 1. Table 3 presents the total control energy over 20 simulation realizations for corresponding  1 .It is noted that for the analyzed model, the control time and the total control energy have meaningful reductions as  1 varies from 0.65 to 0.7.However, when  1 varies from 0.7 to 1.3, and so on, no meaningful reduction of the control time is verified, while the corresponding total control energy has a tremendous increase.
Based on these results, it is possible to come up with the conjecture that for the analyzed model the increase of  1 in (24) significantly reduces the control time and the total control energy up to a limit.Beyond this limit, increase of  1 may not have any further meaningful reduction of the control time while leads to an increase of the total control energy.

Conclusions
The generation of epileptiform spikes in the neural mass models can be derived by increasing the average excitatory synaptic gain and introducing the interactions between populations.The feedback control method based on algebraic estimation is able to control epileptiform spikes in the neural mass models.The influence of the type of the controlled populations, the number of the controlled populations, and the control gain on the control performance is investigated in details for a regular network of three coupled neural populations using a specific control strategy.The results indicate that under the specific control strategy the hyperexcitable populations are necessary choices for being controlled in order to suppress epileptiform spikes.The results also indicate that under the specific control strategy increasing the number of controlled populations leads to the reduction of the control energy, while increasing the control gain does not ensure a meaningful reduction of the control time and the control energy.As a matter of fact, there may exist an optimum value of the control gain that ensures the shortest control time and minimal control energy.
The integral in algebraic estimation of the noisy measurement is evaluated by a trapezoidal numerical integration.The different choices of the integration interval, the sample time, and the number of summation steps lead to different exactness of estimation and different control performance.Therefore, it is important to determine these parameters appropriately for an optimized exactness and control performance.In practice, the integral is able to be approximated by an other method of numerical integration, such as the classical composite Newton-Cotes exact approximation.So, the mathematical theory still offers potential for optimization regarding exactness and control performance, which is only one of the many possible directions of the work.The method of algebraic estimation might approximate not only the signal itself but also the derivative of it from a noisy measurement.There are a variety of fields for which the derivative estimation might be applied.As shown in [30], the derivative estimation is used to construct model-free controller for a shape memory alloy active spring.It should be pointed out that although our algorithm is model based, it may be model-free because the acquisition of []  from algebraic estimation may be independent of system model.No need of a mathematical model is really attractive because a "good" mathematical modeling is difficult if not impossible, especially for the neural system.
In this study, the model parameters are assumed to be homogeneous for populations under consideration.It is not realistic and should be improved in future work.Since the connections of different populations are highly complicated, the influence of topology structure of multiple ( > 3) populations on the control performance is another possible avenue of future work.Validation of the proposed strategy via experimental applications is also one of the possible future directions.

Figure 2 :
Figure 2: Model of multiple coupled neural populations.

Table 1 :
The total control energy over 20 simulation realizations for the model containing one hyperexcitable population.

Table 2 :
The total control energy over 20 simulation realizations for the model containing two hyperexcitable populations. 1 =  2 = 0.86,  3 = 0 and all populations by setting  1 = 0.79,  2 = 0.79,  3 = 0.395.Table2presents the comparison of total control energy over 20 simulation realizations for the above two control schemes.The results indicate that the total control energy decreases with the increase of the number of controlled populations.

Table 3 :
The total control energy over 20 simulation realizations for different values of  1 .