Closed-loop control plays an important role in the treatment of epileptiform spikes by using brain stimulation. In recent years, there have been many analytical methods for determining stimulus protocols and stimulus parameters. However, the analytical method that can start the stimulus protocol when it is needed and stop the stimulus protocol when it is not needed is rather rare. In this work, we design an analytic closed-loop control scheme which can starts control when epileptiform spikes are detected and stops control when no epileptiform spikes are detected. The neural mass model is used to simulate the generation of normal Electroencephalograph signals and epileptiform spikes. The detection of epileptiform spikes is completed via an alarm threshold which is set by using the combination of cross approximate entropy, the Pearson correlation coefficient and the fuzzy theory. If the detection result shows that there are epileptiform spikes in the neural mass model, the fuzzy proportion integration differentiation control works so that the abnormal epileptiform spikes are restored to normal EEG signals, and vice versa. The simulation confirms the effectiveness of the proposed closed-loop control scheme.
Epilepsy is a class of brain disorders caused by abnormal discharges of brain nerve cells. The electrical activity of brain recorded by EEG appears epileptiform spikes with different oscillatory frequencies during seizures [1]. The functional disorders caused by seizures bring physical and mental sufferings and lead to huge family and social burdens [2]. Therefore, the detection and modulation of epileptiform spikes play an important role for the prevention and treatment of epilepsy, as well as for promoting the functional rehabilitation of patients.
In 1954, Penfield and Jasper [3] first described cortical electrical stimulation could reduce seizures. Since then, many therapy methods to alter cortical excitability, such as cortical electrical stimulation of epileptogenic zone [4,5,6,7,8,9,10,11,12], deep brain stimulation [13,14,15,16,17], transcranial repetitive magnetic stimulation [18,19,20], vagus nerve stimulation [21,22,23,24] and trigeminal nerve stimulation [25, 26], have been available to normalize disrupted forms of brain activity for patients with epilepsy. Brain stimulation seems to be gradually replacing surgical surgery as the most promising means of treatment of epilepsy. However, both the electrical stimulation and the magnetic stimulation, where the selection and adjustment of the stimulus parameters such as frequency, amplitude and pulse width, are largely dependent on the actual experience of the designer or surgeon. The therapeutic effect of brain stimulation is far from optimal. The closed-loop stimulation constructed by feedback of EEG signals can detect the abnormal electrical signals of lesion sources in real time, and adjust adaptively the stimulus parameters according to the specific situation of detection, so that the therapeutic effect is optimized, and the side effect, the degree of damage and the energy loss are reduced [27, 28]. In our previous work, we constructed a closed-loop control strategy based on the neural mass model and the unscented Kalman filter algorithm, and confirmed the effectiveness of the strategy in the suppression of epileptiform spikes [29]. Meanwhile, we combined the closed-loop control strategy with the fuzzy theory to further optimize the control performance [30]. In view of the complexity and high nonlinearity of EEG signals, we designed a fuzzy PID control strategy based on the neural mass model to make the output track an expected waveform [31]. The feedback linearization control was applied to suppress seizures based on the neural mass model in [32].
Such closed-loop control frameworks have a common feature, that is, the control action is always imposed on the controlled object regardless of seizures or not. The actual situation is that epilepsy patients are not always in the state of seizures. It is not necessary to exert control when there is no seizures. Otherwise it may cause unnecessary energy consumption and other problems, and in recent years, the use of entropy to classify EEG signals has become a more and more important analysis method. In view of this, we design a new closed-loop control framework based on the neural mass model which combines cross approximate entropy, the Pearson correlation coefficient, the PID control algorithm and the fuzzy theory together to detect epileptiform spikes and impose the control action on demand according to the detection results. Specifically, the control action is imposed on the controlled object and epileptiform spikes are modulated to normal EEG signals when the detection results confirm that there exhibit epileptiform spikes in the controlled object. On the contrary, the control action stops working when the detection results confirm that there are no epileptiform spikes in the controlled object.
Neural mass model. Fig. 1 presents the structure and block diagram of the neural mass model with multiple coupled populations. As shown in Fig. 1 A, the individual neural mass model mainly consists of two subsets. One denotes pyramidal cells, and the other denotes local interneurons (i.e. other non-pyramidal cells, stellate or basket cells). Pyramidal cells receive excitatory and inhibitory feedback from interneurons and excitatory input from neighboring or more distant populations. Local interneurons receive excitatory input only from pyramidal cells. The interaction between pyramidal cells and local interneurons produces oscillations. Two functions named "pulse-to-wave function" and "wave-to-pulse function" by freeman [33, 34] are used to characterize each subset. The former is a linear transfer function that transforms presynaptic information (average density of afferent action potentials) into postsynaptic information (average excitatory or inhibitory postsynaptic membrane potential). The latter is a static non-linear function that associates the average level of membrane potential of the neurons with an average pulse density of potentials fired by these neurons. The linear transfer function can be represented by a second order low pass filter with the form:
Fig. 1.(A) The structure of the neural mass model. (B) Block diagram of the neural mass model.
$ H(s) = G/(s + g)^2, $
where $ s $ is the Laplace variable, $ G = A $ and $ g = a $ for the excitatory case, $ G = B $ and $ g = b $ for the inhibitory case, $ A $ and $ B $ stand for the average excitatory and inhibitory gain, $ a $ and $ b $ are the membrane average time constant and dendritic tree average time delays. The parameters $ A $, $ B $, $ a $ and $ b $ can modulate the sensitivity of excitatory and inhibitory synapses. The impulse response of the linear transfer function is given by:
for the excitatory case and:
for the inhibitory case. The static non-linear function is represented by a sigmoid function:
where $ e _{0} $, $ v _{0} $ and $ r $ represent the maximum firing rate, the postsynaptic potential corresponding to $ e _{0} $ and the steepness of the sigmoid function, respectively. A gain constant $ {K }_{lj} $ is applied to define the degree of coupling between population $ l $ and population $ j $, where $ l = 1, 2, \cdots, M$, j = 1,2,$\cdots$,M, $ {l\neq j} $, $ M $ is the total number of considered populations. The corresponding block diagram of population $ l $ ( l = 1,2,$\dots$,M) in the neural mass model is accordingly derived, as shown in Fig. 1 B. The model can be considered as a feedback system driven by a noise input $ p ^1 (t) $ that globally denotes the average density of afferent action potentials from neighbouring or distant populations. Here p$ ^1 (t)$ can be any arbitrary function including white noise. The delay associate with connections from population $ l $ is characterized by a filter with an impulse response $ h_d (t) $, $ h_d (t) $ be chosen as being similar to $ h_e (t) $:
where $ a_d $ represents the average time delay on efferent connection from a population. The variables $ y^l_0 (t) $, $ y^l_1 (t) $ and $ y^l_2 (t) $ denote the feedback of the pyramidal cells to the excitatory postsynaptic membrane potential of the interneurons, the average excitatory and inhibitory postsynaptic membrane potential, respectively. The substraction of $ y^{l}_2 (t)$ from $ y^l_1 (t) $, i. e., $ y ^{l}_1 (t) - y ^{l}_2 (t) $, is to be taken as representative of the cortical EEG. Let $ y^l (t) = y^l _1 (t)- y^l_2 (t) $. The variable $ y^l_6 (t) $ denotes the output of $ h_d (t) $ which establishes the connection from population $ l $ to other populations. The interaction between the pyramidal cells and interneurons is characterized by four connectivity constants $ C_{1} - C_{4} $, which account for the average number of synaptic contacts.
According to (1), (2), (4) and Fig. 1 B, we can deduce that each function ($ h_e (t) $, $ h_i (t) $ or $ h_d (t) $) introduces a two-order differential equation
Letting:
the formulas in (5) can be transformed into:
In conclusion, the neural mass model given in Fig. 1B can be represented as the following set of eight differential equations:
All parameters in the neural mass model are set on a physiological basis, and their default values are:
Detailed information on how these values were obtained is given in [35].
Cross approximate entropy. Consider two time series $ \{u(1), u(2)$, $ \dots, u(N) \}$ and $ \{ v(1), v(2), \dots, v(N) \}$. The algorithm of cross approximate entropy is given as follows.Firstly, two $ m- $dimensional vectors $ U(i )$ and $ V (j)$ are constructed by using the components of the two time series given above
Secondly, the distance between $ U(i) $ and $ V(j) $ is defined as:
Set a threshold value $ r > 0 $. Calculate the number of $ d [U(i), V (i)] < r $ as $ C $, the ratio of the $ C $ to the total number of vectors is given as $ C^m_i (r) $:
Thirdly, take the logarithm of $ C^m_i (r) $ and denote the average of ln $ C^m_i (r) $ for all $ i $ as $ \phi m(r) $:
Similarly, $ \phi m+1 (r) $ can be obtained by repeating above steps. Finally, The cross approximate entropy of the two time sequences is defined as:
$C-ApEn (m, r) = \phi m (r) - \phi m+1 (r)$
Pearson product-moment correlation coefficient. For two given time series $ \{ u(1), u(2), \dots, u(N) \}$ and $ \{ v(1), v(2), \dots, v(N) \}$, the mathematical expression of pearson product-moment correlation coefficient (PPMCC) is defined as:
where $ u $ and $ v $ are the average values of the time series $ \{u(1), u(2)$, $ \dots , u(N)\}$ and $\{v(1), v(2)$, $ \dots , v(N)\}$, respectively. The value of $ P_{u,v} $ lies between $ -1 $ and $ +1 $, i. e., $ |P_{u,v}| \leq 1 $. If $ |P_{u,v}| $ is greater, the correlation between $\{u(1), u(2), \dots , u(N)\}$ and $ \{v(1), v(2), \dots , v(N)\}$ is stronger.
In general, long-term observation of people with epilepsy is required as seizures are unpredictable random episodes. The dynamic characteristics of EEGs can provide important information for diagnosis and the monitoring of epilepsy. However, the method of visual inspection is time-consuming, has low efficiency, and lacks standards. Automatic detection of seizures during long-term EEG records is much preferred. In previous work, traditional PID control has been combined with FIS to form a fuzzy PID control algorithm to modulate epileptiform spikes generated by the neural mass model. It achieves this by taking into account the nonlinearity and high complexity of the model [31]. On the basis of the closed-loop framework given in [31], a new closed-loop control scheme was designed for the neural mass model. This enables on demand modulation of epileptiform spikes according to the epileptiform spikes detected.
Fig. 2 presents a diagram of the closed-loop control scheme designed for the neural mass model. The whole diagram can be divided into three parts. The first part is the controlled object, i. e., the neural mass model. The symbols $ y ( t ), y ^{\ast} (t) $ and $ e(t) = y^{\ast} (t) - y ( t) $ denote the column vector constructed from the output of three populations, the expected output which is viewed as the desired modulation, and the error signal which represents the degree of deviation of $ y (t)$ from $ y^{\ast} (t) $, respectively. The second part is the fuzzy PID control algorithm imposed on the neural mass model. The symbol $ u(t) $ denotes the control law and has the form of the fuzzy PID control:
Fig. 2.Diagram of the closed-loop control scheme.
where $ k_p $ , $ k_i $ and $ k_d $ are the proportion coefficient, the integral coefficient and the differential coefficient, respectively. These coefficients are derived from online adaption system, the principle of which is the same as that given in [31]. The third part, marked by the dashed block, is used to detect whether epileptiform spikes exhibited and decides whether the switch is closed. This part is mainly composed of PPMCC method, Cross approximation entropy, FIS, multiplier, comparator, and the alarm threshold. The inputs of the blocks "Cross approximate entropy" and "PPMCC method" are $ y (t)$ and $ y ^{\ast} (t) $, more specifically, the output of the hyperexcitable population 1 and the component $ y_1^\ast (t) $ of the expected output. The outputs of the blocks "Cross approximate entropy" and "PPMCC method" are $ p $ and C - ApEn respectively. The block "FIS" is introduced to distinguish the normal EEGs and epileptiform spikes better, The inputs of FIS are $ p $ and C - ApEn and the output of FIS is $ S $. The fuzzy inference system is mainly composed of fuzzification, a fuzzy rule base, a fuzzy reasoning method, and defuzzification.
Although there are no restrictions on the form of the input membership functions, the triangular membership function is used in the premise mainly because it has the characteristics of less computation, easy implementation and good performance. A triangular membership function is specified by three parameters ($ a $, $ b $, $ c $, ) and has the following form:
where $ {a \leq b\leq c} $, $ \mu $ is the membership grade, $ b $ represents the center of the membership functions, and $ a $ and $ c $ determine the endpoint. The Mamdani minimax reasoning method is used for the fuzzy logic reasoning. The membership functions are normalized using five membership functions NB, NM, Z, PM, PB, as show in Fig. 3. Here "NB" expresses the negative big, "NM" expresses the negative medium, "z" expresses the zero, "PM" expresses the positive medium, "PB" expresses the positive big. To contain the value accurately, z-shaped and s-shaped membership functions are used as the boundary curve of $ P $ for PB and NB. The generation methods of fuzzy control rules are based on the experience and knowledge of experts, operator's mode of operation and learning algorithms. In the decision-making process, the posteriori fuzzy inferences are determined by the twenty rules given in Table 1. These rules mainly control the degree of influence of $ S $ on $ C $ - ApEn. Finally the gravity center method is used for defuzzification and derivation of $ S $. The multiplier of $ S $ and $ C $ - ApEn:
is used for comparisons with the alarm threshold $ \theta $. Once the comparison of results shows the presence of epileptiform spikes, the switch is closed and the fuzzy PID control is imposed on the neural mass model such that epileptiform spikes are modulated into normal EEG signals. Otherwise, the switch is open and fuzzy PID control is idle.
Fig. 3.Input membership function.
|
An important aspect of the detection of epileptiform spikes is the setting of an appropriate alarm threshold. This is achieved by using the method given in Fig. 2. It has been shown that the main cause of seizures is increased excitability of a a brain region [1]. In the neural mass model, the ratio A/B controls the degree of excitability within a given population. If $ A $ is increased from the standard value while other parameters are maintained at their standard values, the neural mass model produces epiletiform spikes. Thus, the alarm threshold $ \theta $ can be set based on these features. Consider the neural mass model coupled to three neural populations as shown in Fig. 4. The value of $ A $ of population 1 varies from 3.2mV to 4mV and the remaining parameters are set at the standard values. The coupling coefficients $ K_{12} $, $ K_{23} $ and $ K_{31} $ are set at 100. The function $ p ^1 (t) $ is chosen as a Gaussian white noise, mean value 101, and standard deviation 35.
Fig. 4.Neural mass model coupled by three neural populations.
Fig. 5.Relationships between $ \tilde\theta $ - and $ A $.
Fig. 6.Output of the neural mass model.
Fig. 7.Relationships between $ C $ $ \mathrm{-} $ ApEn and $ A $.
The corresponding time series of $ y_{1} (t) $ and $ y _1^{\ast} (t) $ is chosen as $ u(i) $ and $ v(i)$, and time step 0.001s. During calculation of C - ApEn, the embedding dimension $ m $ is selected as two and the threshold value $ r $ as 0.2cov[$ u(i), v(i)$]. The setting of the values of $ m $ and $ r $ gives reasonable statistical characteristics for judging complex signals [36]. Fig. 5 gives the relationships between $ \tilde\theta $ and $ A $. The value of $ \tilde\theta $ increases immediately from a small value to a relatively large value when $ A $ equals 3.33mV. The alarm threshold $ \theta $ is set to 0.1 by taking into account the relationships between both the dynamics and the value of $ A $ and the relationships between $ \tilde\theta $ and $ A $. In other words, epileptiform spikes exist in the controlled object if $ \tilde\theta \geq $ 0.1. Fig. 6 gives the reference output of the neural mass model. It exhibits no epileptiform spikes, sporadic random epileptiform spikes, and sustained spike discharges when the value of $ A $ of population 1 is set to 3.25mV, 3.34mV, or 3.6mV, respectively. These dynamics reflect normal activities that resemble those of real EEG signals and real EEG signals during the propagation of temporal lobe seizures [1]. For comparison, the relationship between C - ApEn and $ A $ is shown in Fig. 7. It was observed that C - ApEn gradually increases with fluctuations when $ A $ increases. The value of C - ApEn is not very different between adjacent pairs of points. It is difficult to determine an appropriate threshold according to Fig. 7. Thus, the combination of the PPMCC and the FIS with cross approximate entropy well implements the detection of epileptiform spikes in the controlled object. To further confirm the validity of the given threshold selection method, more simulations are provided. In Fig. 4, the parameter settings are the same as for Fig. 5, with the exception that the inhibitory gain $ B $ of all populations is set to 20mV and 24mV, respectively. Fig. 8 presents the relationship between $\tilde \theta $ and $ A $, where the green curve and the pink curve are the results for $ B $ = 20mV and $ B $ = 24mV, respectively. Comparing the results in Fig. 8, it was concluded that the relationships between $ \tilde\theta $ and $ A $ under different settings of $ B $ have similarities, that is, $ \tilde\theta $ suddenly increases from a small value to a relatively large value for certain values of $ A $. Thus, it is reasonable to set the alarm threshold $ \theta $ to 0.1.
Fig. 8.Relationships between $ \tilde\theta $ and $ A $.
Fig. 9.Dynamics of the neural mass model without and with the action of control.
The neural mass model given in Fig. 4 is considered to verify the effectiveness of the proposed closed-loop control scheme. To simulate the situation that epileptic patients are normal for most of the time and in a seizure state for only short durations, the value of $ A $ of population 1 is set randomly from 3.2mV to 3.4mV with high probability and set greater than 3.4mV with low probability. The other parameters are set to their default values, and the coupling coefficients $ K_{12} $, $ K_{23} $, and $ K_{31} $ are set to 100. The dynamics of the neural mass model under these parameter settings are shown in the left column of Fig. 9. It is observed that the controlled object exhibited epileptiform spikes at 58s that continue thereafter. The moving window method is used to calculate $ \tilde\theta $ in real time and the window moves 1s each time. The blue curve at the bottom of Fig. 9 denotes the change of $ \tilde\theta $ with time. The value of $ \tilde\theta $ increases suddenly from small values to a relatively large value at 58s, thus indicating initiation of epileptiform spiking in the controlled object at that time. This detection result is consistent with the actual case of the controlled object. Thus, the detection method of epileptiform spikes is also efficient for the case of the neural mass model actually considered. As the closed-loop control scheme shows, the switch is closed and fuzzy PID control is imposed on the neural mass model once epileptiform spikes are detected. The right column of Fig. 9 is the modulated result after activating the detection and control components. The output with epileptiform spikes is converted into a normal EEG signal within about 3s. The change of $ \tilde\theta $ with time at the bottom of Fig. 9 (green curve) confirms this point.
Fig. 10.Network model coupled to eight neural populations.
Fig. 11.Dynamics of the neural mass model with and without control action.
Further simulations are provided to verify the effectiveness of the closed-loop control scheme. Consider a network model coupled by eight neural populations, as shown in Fig. 10. The value of $ A $ of population 1 is randomly set with either high probability between 3.2mV and 3.25mV or with a low probability of being set greater than 3.4mV. The value of $ B $ of all populations is set at 24mV. The remaining parameters are set at the standard values and the coupling coefficients are set to 30. The dynamics of the network model under these parameter settings are shown in the left column of Fig. 11. It was observed that the controlled object exhibits epileptiform spikes at 38s that continue after that time. The time window at the bottom of Fig. 11 shows that $ \tilde\theta $ suddenly jumps from a small value to a relatively large value at 38s. Thus, the detection of epileptiform spikes is in agreement with that of the controlled object, and the detection method for epileptiform spikes is efficient for the network model. The right column of Fig. 11 presents the resulting dynamics of the network model following activation of closed-loop control. Epileptiform spikes disappear within a few seconds and the change of $\tilde \theta $ with time at the bottom of Fig. 11 (green curve) confirms this.
In summary, fuzzy PID control has advantages of speed and stability in the modulation of epileptiform spikes produced by the neural mass model. Meanwhile, the entire closed-loop control scheme, including the detection of epileptiform spikes and activation of fuzzy PID control, is validated for the automatic detection and modulation of epileptic EEG signals produced by the neural mass model.
Fig. 12.Dynamics of the neural mass model with and without control action.
Fig. 13.Variation of control energy with time and the total control energy for the proposed closed-loop control scheme given in Fig. 3.
Fig. 14.Variation of the control energy with time and the total control energy for the conventional fuzzy PID control.
An important consideration in the modulation of epileptiform spikes is the control energy required to achieve the modulation process. The proposed closed-loop control scheme has advantages in the control of energy because of the presence of epileptiform spike detection. In what follows, this point is demonstrated by providing simulation data for the neural mass model, given in Fig. 4. The simulation time was 40s and evenly divided into four periods. The value of $ A $ for population 1 during each period was set to 3.25mV, 3.26mV, 3.44mV, and 3.29mV. Other parameters are set to their default values and the coupling coefficients $ K_{12} $, $ K_{23} $, and $ K_{31} $ are set to 100. The dynamics of the neural mass model under these parameter settings exhibits epileptiform spikes in $ 20 s\sim 30s $ and normal EEG signals in the remaining time, as shown in the left column of Fig. 12. The features of the dynamics can be deduced by analysis of results and used to determine the alarm threshold. To verify the advantages of the control energy required to implement the epileptiform modulation, a conventional fuzzy PID control that lacked the detector and the proposed closed-loop control scheme given in Fig. 2 are used respectively for modulation of abnormal EEG signals. The control energy is defined as $ u^T (t) u (t)$. The right column of Fig. 12 presents the result of using the proposed closed-loop control scheme given in Fig. 2, where the blue curve denotes the desired normal EEG signal and the red curve denotes the output of the neural mass model the control is applied to. It shows that the proposed closed-loop control scheme was successful in achieving the appropriate modulation. The output curve of the neural mass model with the imposition of conventional fuzzy PID control without the detection component is similar to that given in Fig. 12 and thus is omitted. Fig. 13 and Fig. 14 illustrate the variation of the control energy with time and the total control energy required to realize the required modulation for the conventional fuzzy PID control and the proposed closed-loop control scheme. Fig. 14 shows that the control energy is 0mV during the first, second, and fourth intervals during which the controlled object exhibits normal EEG signals, and is 192.9 $ \times$ 103mV2 during the third interval when the controlled object exhibited epileptiform spikes. However, it can be found from Fig. 13 that the control energy is always generated whether or not there are epileptiform spikes in the controlled object. The total control energy needed is 235.5 $ \times $103mV2. Thus, the proposed closed-loop control scheme has advantages when compared to the control energy otherwise required for successful modulation.
A new closed-loop framework has been proposed here that enables a neural mass model to modulate abnormal epileptiform spikes into normal EEG signals. The most obvious feature of the proposed framework is its ability to control on demand according to the detection of epileptiform spikes. In other words, the controller acts on the neural mass model when seizure is detected and stops working on the neural mass model in the absence of seizures. Thus, the advantages of the proposed framework have been demonstrated for the control energy required to achieve appropriate modulation.
In this work, the detection of epileptiform spikes and the fuzzy PID control are combined to achieve modulation in an EEG generation model initially proposed by Lope Da Silva [33] and extended by Wendling et al. [1]. Recently, more sophisticated models that can simulate EEGs with different rhythms have been proposed such as the model of impaired gabaergic dendritic inhibition [37], realistically coupled neural mass models [38], the neural mass model for investigation of event-related potentials in the frequency domain [39], nonlinear conductance-based neural population models [40], a model to fit the impulse response in three cortical regions of interest recorded during a series of TMS/EEG experiments [41], and a lumped computational model of the thalamo-cortico-thalamic circuitry [42]. It is planned that some of these sophisticated models will be incorporated to study the detection and modulation of spikes. The alarm threshold is obtained by using cross approximation entropy and an adjustment factor derived from a fuzzy inference system. The current work can also be extended to another closed-loop scheme presented in earlier work [29], where an unscented Kalman filter control is used to suppress the epileptiform spikes of a neural mass model. It is of interest and value to extend the current framework to networks of neural populations for the modulation of brain dynamics. It is hoped this work can further assist in the treatment of epilepsy and other brain diseases.
This work was supported by National Natural Science Foundation of China (61473245, 61004050), and by Natural Science Foundation of Hebei (F2017203218).
All authors declare no conflicts of interest.