A Method to Obtain Parameters of One-Column Jansen–Rit Model Using Genetic Algorithm and Spectral Characteristics

: In this paper, a method of obtaining parameters of one-column Jansen–Rit model was proposed. Methods present in literature are focused on obtaining parameters in an on-line manner, producing a set of parameters for every point in time. The method described in this paper can provide one set of parameters for a whole, arbitrarily long signal. The procedure consists of obtaining speciﬁc frequency features, then minimizing mean square error of those features between the measured signal and the modeled signal, using genetic algorithm. This method produces an 8-element vector, which can be treated as an EEG signal feature vector speciﬁc for a person. The parameters which were being obtained are maximum postsynaptic potential amplitude, maximum inhibitory potential amplitude, ratio of the number of connections between particular neuron populations, the shape of a nonlinear function transforming the average membrane potential into the ﬁring rate and the input noise range. The method shows high reproducibility (intraclass correlation coefﬁcient for particular parameters ranging from 0.676 to 0.978) and accuracy (ranging from 0.662 to 0.863). It was additionally veriﬁed using EEG signal obtained for a single participant. This signal was measured using Emotiv EPOC+ NeuroHeadset.


Introduction
The electroencephalography (EEG) is a noninvasive brain examination, in which brain waves are measured on a skin surface of the subject's head. EEG is most often used by neurologists to differentiate functional from organic brain diseases, to diagnose sleep disorders, headaches, dizziness, to monitor brain activity during heart operations. EEG offers high temporal resolution which is not possible with magnetic resonance imaging (MRI).
Since the EEG's history goes back to 1924 [1], the EEG signal analysis field is profoundly studied [2][3][4]. On the other side, there is whole field devoted to EEG signal modeling, which can be divided into three main categories [5]. Firstly, there are microscopic (detailed) models, focused on modeling individual neurons and their connections. Secondly, mesoscopic (reduced) models, in which neurons are modeled in less complex manner. Finally, macroscopic models, which simulate whole neuron networks without simulating individual neurons. In this paper, the third category is considered.
First macroscopic model, considering only excitatory neural interactions, was introduced in 1956 by Beurle [6]. In 1972, Wilson and Cowan proposed a model with additional inhibitory interactions [7]. Amari, in 1975, proposed a spatial-temporal model [8], in opposition to only time-variant predecessors. Despite the simplicity of those models, they constitute the foundation of all the later ones [5], like Lopes da Silva's [9], Jansen's [10], or Wendling's [11].
In this paper, the one-column Jansen-Rit (J-R) model was analyzed. The mathematical single-column model was developed by many researchers including Lopes Da Silva, Jansen, and Rit [10]. It is widely present in literature [10,[12][13][14][15][16][17][18][19][20][21][22]. In [10], the authors deeply analyzed the model, taking into account not only alpha waves, but also epileptic activity and bifurcations, with the emphasis put on the influence of the p parameter (so the representation of the activity exterior to the column). In [12], the authors focused on comparison between the one column J-R model, double column J-R model, and Wendling's model. However, parameters of the models were adjusted to represent specific brain activities, without comparison to real signals. Authors of [14] extended the model and analyzed both real and simulated signals, with the focus on transition from sine to spike activities. Nonetheless, fitting the model to real signal was absent. In [16,17], the J-R model was extended, however without explicit, numerical comparison to real signals. The one-column J-R model has been also used under stochastic and periodic driving, to study ictogenesis (progression from the interictal state to the seizure) [19]. In [20], the authors studied the dynamical interactions of cortical areas. In [21], the authors used Jansen-Rit model to develop an in silico platform to examine the effects of neighboring regions electrical stimulation on a another brain region. The authors of [22] analyzed the dynamics of a revised Jansen-Rit model in the stop-signal task.
In 2012, the authors of [23] presented a method for estimating parameters in the J-R model using state observer. In the next year the method was implemented to estimate parameters of both modeled and real signals (which were obtained via subdural grid electrode) [24]. In 2016, the authors of [25] applied unscented Kalman filter to estimate the J-R model parameters in real time, while monitoring propofol-based anesthetic state (via EEG). Obtained parameters were used to classify anesthetic brain state.
In the mentioned works, the parameters were being obtained in a continuous way; i.e., they constituted a parameter-time matrix, with a set of parameter values for each specific moment. This paper presents a different approach to obtaining J-R model parameters, which is an off-line method. Obtained parameters, like in previous works, could be used for signal classification. The parameters themselves, however, are obtained once for full duration of a signal. That means that they could hypothetically constitute EEG signal features for specific individual, enabling off-line classification of obtained signals.
To obtain the parameters for the whole signal, the genetic algorithm was used. Minimized function was defined as the mean square error (of the modeled and real signal) of the power spectrum of some specific frequencies. The specific objective of this research was to test if minimizing function defined in this way would be suitable, as it proved to be. Empirical comparison shows that the generated signals show close similarity to measured signals in the spectrum content.

Modeled Signal
The conducted research included modeling of the neural structure by the Jansen-Rit single column model [10]. It includes a population of pyramidal structures receiving stimulatory and inhibitory impulses from interneurons and an input fed with noise simulating the impact of other neural structures. Pyramidal cells of the cerebral cortex, due to their characteristic location in the cerebral cortex structure, are considered to be the main source of the electroencelographic signal from the human brain. The Jansen-Rit model's block diagram is shown in Figure 1.
The P 1 , P 2 , and P 3 blocks can be presented using the equations:  (1), P3 is the function (hi) presented in Equation (2). The S1, S2, and S3 (sigm) blocks are the sigmoidal functions presented in Equation (5), p is noise, while the C1-C4 are the coefficients of connections between particular neural structures.
The P1, P2, and P3 blocks can be presented using the equations: where: A and B are the maximal post-synaptic potential amplitudes, a and b are membrane time-constants. Equation (1) can also be presented in the form of two first order differential equations: Equation (2) also can be expressed this way, by changing "A" and "a" terms to "B" and "b". The model can be therefore expressed using six first-order equations [10]: where sigm is the sigmoidal function, responsible for the neural structure activation rate, dependent on the potential (graphical representation of this function was presented in Figure 2): p(t) is the signal delivered to the model, made of white noise, originally in the range of 120-320 pulses per second, C1, C2, C3, and C4 are numbers of synapses established between particular neuron populations, bounded by equation: Figure 1. The one-column Jansen-Rit model. The P 1 and P 2 blocks are the functions (he) presented in Equation (1), P 3 is the function (hi) presented in Equation (2). The S 1 , S 2 , and S 3 (sigm) blocks are the sigmoidal functions presented in Equation (5), p is noise, while the C 1 -C 4 are the coefficients of connections between particular neural structures.
Equation (1) can also be presented in the form of two first order differential equations: .
In the sigmoidal function the parameters e 0 and v 0 are directly responsible for the activation function's shape, i.e., the inflection point. The parameter r is responsible for Appl. Sci. 2021, 11, 677 4 of 12 steepness (slope). These parameters were being modified in a course of this study, since the postsynaptic potential transforms into a firing rate not exactly in the same way from a population to a population [10,26].
The model was simulated with frequency of 1 kHz, using mid-point numerical method of solving differential equations. Output of the model, or the modeled EEG signal, is denoted as y and can be obtained by simple subtraction: In the sigmoidal function the parameters e0 and v0 are directly responsible for the activation function's shape, i.e., the inflection point. The parameter r is responsible for steepness (slope). These parameters were being modified in a course of this study, since the postsynaptic potential transforms into a firing rate not exactly in the same way from a population to a population [10,26].
The model was simulated with frequency of 1 kHz, using mid-point numerical method of solving differential equations.

Model Parameter Obtaining Process
Due to the small variability in the time constant parameters of the membrane and delays in the dendritic tree (a and b, [10]), it was assumed that the model can feature changes in eight parameters, i.e., A, B, C, e0, v0, r, and p's lower limit and range. These parameters mostly affect the form of the signal. Due to the relatively high number of arguments in the minimized function, the genetic algorithm (GA) was used for the minimization process. Along with other heuristics, it is often used for parameter estimation [27].
Preliminary study was done comparing GA to different common heuristics, i.e., simulated annealing (SA), particle swarm (PS), and surrogate optimization (SO). In the task of J-R model parameter obtaining, the GA and the PS proved to be more accurate than the rest of the compared algorithms. Despite the longest time to compute, genetic algorithm was chosen because of its lowest value of found optimum and vast tuning capabilities. Specific values of the computation time and the found optimum are presented in Table 1. Note, that both, time and optimal value, were normalized to the highest value.

Model Parameter Obtaining Process
Due to the small variability in the time constant parameters of the membrane and delays in the dendritic tree (a and b, [10]), it was assumed that the model can feature changes in eight parameters, i.e., A, B, C, e 0 , v 0 , r, and p's lower limit and range. These parameters mostly affect the form of the signal. Due to the relatively high number of arguments in the minimized function, the genetic algorithm (GA) was used for the minimization process. Along with other heuristics, it is often used for parameter estimation [27].
Preliminary study was done comparing GA to different common heuristics, i.e., simulated annealing (SA), particle swarm (PS), and surrogate optimization (SO). In the task of J-R model parameter obtaining, the GA and the PS proved to be more accurate than the rest of the compared algorithms. Despite the longest time to compute, genetic algorithm was chosen because of its lowest value of found optimum and vast tuning capabilities. Specific values of the computation time and the found optimum are presented in Table 1. Note, that both, time and optimal value, were normalized to the highest value. Function minimized in this comparison, and in the rest of this paper, is described below.

Function Minimized in This Study
In order to match the EEG signal generated by the model to the reference EEG signal, it was necessary to develop an error function that would be minimized by the genetic algorithm.
The frequency analysis was used in the function to limit the impact of the differences in length between the modeled and real signals. The value of the minimized function was calculated as the mean square error (of the modeled and real signal) from the power spectrum obtained via discrete Fourier transform (DFT) in range 2-18 Hz. Note that the dimensions of this product may change depending on the sampling frequency and the length of the signal.
Function defined in that way is dimensionless. Its visualization is given in Figure 3. Function minimized in this comparison, and in the rest of this paper, is described below.

Function Minimized in This Study
In order to match the EEG signal generated by the model to the reference EEG signal, it was necessary to develop an error function that would be minimized by the genetic algorithm.
The frequency analysis was used in the function to limit the impact of the differences in length between the modeled and real signals. The value of the minimized function was calculated as the mean square error (of the modeled and real signal) from the power spectrum obtained via discrete Fourier transform (DFT) in range 2-18 Hz. Note that the dimensions of this product may change depending on the sampling frequency and the length of the signal.
Function defined in that way is dimensionless. Its visualization is given in Figure 3.

The Genetic Algorithm
To minimize the described function (MSE of the power spectrum), the Global Optimization Toolbox of the Matlab environment was utilized [28].
To ensure robustness and avoid possible local minima problems, the algorithm was conducted with the population size of 256. The parameter obtaining process stopped after 150 populations, when the performance stagnated. Selection of parents for the consecutive generations was conducted using the stochastic universal sampling. In the reproduction step, 13 (5% of the population size) best individuals were guaranteed to survive to the next generation. The remaining 243 individuals were obtained in 80% by crossover and in 20% by Gaussian mutation.
The algorithm sought the function values in the ranges presented in Table 2, among real numbers.

The Genetic Algorithm
To minimize the described function (MSE of the power spectrum), the Global Optimization Toolbox of the Matlab environment was utilized [28].
To ensure robustness and avoid possible local minima problems, the algorithm was conducted with the population size of 256. The parameter obtaining process stopped after 150 populations, when the performance stagnated. Selection of parents for the consecutive generations was conducted using the stochastic universal sampling. In the reproduction step, 13 (5% of the population size) best individuals were guaranteed to survive to the next generation. The remaining 243 individuals were obtained in 80% by crossover and in 20% by Gaussian mutation.
The algorithm sought the function values in the ranges presented in Table 2, among real numbers.

Computation Complexity
Since genetic algorithm is a stochastic approach to solve a problem, it is, by definition, impossible to exactly determine its computation complexity. However, with restraints imposed on some of its parameters, the maximal computation complexity (in terms of time) can be determined.
In the proposed method, the main restraint is the number of populations. After reaching this number, the algorithm stops, and the best score obtained by it is given as a result. Preliminary research pointed out that for the defined fitness function (MSE of the power spectra), the algorithm convergences in less than 150 iterations. Therefore, this number was chosen as a maximal number of populations.
Knowing that every population is composed of a constant number of 256 individuals, it can be determined that for every genetic algorithm execution there will be maximally 38,400 (150 populations, each made from 256 individuals) fitness function evaluations. To determine the time needed for single fitness function evaluation, the function was executed for 100,000 times with different, random inputs. The mean time of those executions was used to determine maximal computation complexity of a proposed method.
The genetic algorithm itself can be sped up by various means; for example, saving value for each individual and just loading the value instead of evaluating it, or imposing additional restraints, like, for example, maximum number of stall generations after which the algorithm stops. The value determined in this study is the maximum time needed for the complete GA execution.
The study of time needed for the minimized function execution was conducted on the computer with Intel Core i7-3740QM CPU (four 2.7 GHz cores, although function was evaluated without parallelization, i.e., using only one core), Crucial MX500 SSD drive (560 MB/s read speed), and 16 GB RAM (4 × 4 GB DDR3, 1600 MHz), in the Matlab R2020a environment.

Repeatability of the Method
To determine the repeatability of the proposed method, the intraclass correlation coefficient (ICC) was used. To obtain the ICC, k measurements should be conducted on each of n subjects. In this study, the parameter obtaining process should be interpreted as a measurement. According to [29], 100 subjects with 10 measurements per subject should be more than sufficient for the accurate ICC determination.
The ICC can be defined in number of different ways, depending on the use. Here, the coefficient was defined for "two-way mixed effects, absolute agreement model" (in literature known as "ICC(A,k)" [30]), as suggested in [31]. It is defined as follows: where: MSBS is the mean square between modeled signals, MSE is the mean square error, MSBM is the mean square between measurements (parameter obtaining processes), n is the number of the measured signals.

Accuracy of the Method
As mentioned before, to the best knowledge of the authors, in the literature, there are no deterministic ways to obtain the J-R model parameters for the real (measured) EEG signals. Therefore, it is impossible to compare parameters obtained via the proposed method to the "real" parameters of measured signals, as those are unobtainable. The accuracy of the method can be determined, however, using signals generated by the J-R model itself. This way, the original parameters of the signal are known and the difference between them and the parameters obtained by the proposed method could be obtained.
Accuracy was defined as one minus the difference between the obtained parameter's value and the real parameter's value divided by the range of possible parameter's values (from Table 2). It was visualized in the right part of Figure 4. Accuracy defined that way is an 8-element vector with values ranging from 0 (when the given parameter is as wrong as it could be) to 1 (when the parameter's value is identical to the reference).
is an 8-element vector with values ranging from 0 (when the given parameter is as wrong as it could be) to 1 (when the parameter's value is identical to the reference).

Reference Signals Generation Process
Both, the accuracy and the repeatability studies of the proposed method, were conducted on the dataset generated by the J-R model.
First, 100 signals were generated with random parameter values, from ranges listed in Table 2. Then, to somehow simulate real measurement, additional white Gaussian noise was added to each signal. Signal to noise ratio was equal to 0 dB for every signal, which means that the power of added noise was equal to the power of the original signal. Next, for each noised signal, the parameters were obtained 10 times using the proposed method.
This process is graphically presented in Figure 4. For clarity, only one signal generation was visualized. Each reference signal was generated with 1 kHz sampling frequency using a midpoint numerical method of solving differential equations and lasted for 20 s.

Real, Measured Signal
The proposed method was additionally tested on a real signal. To measure it, a 14channel Emotiv EPOC+ NeuroHeadset electroencephalograph was used. Emotiv EPOC+ NeuroHeadset is designed for scalable and contextual human brain research. Arrangement of specific electrodes utilize 10-20 standard. The signal was measured between the O2 electrode and the Common Mode Sense (CMS) electrode [32], which was chosen because alpha waves have biggest amplitude in occipital lobe, especially in the nondominant hemisphere [33]. Signal was sampled with 128 Hz frequency. Obtaining electrode resistance not exceeding 10 k Ohm at the start of the measurement is an essential condition Figure 4. Visualization of the accuracy study. Note, that this process was conducted 100 times, with different reference parameters.

Reference Signals Generation Process
Both, the accuracy and the repeatability studies of the proposed method, were conducted on the dataset generated by the J-R model.
First, 100 signals were generated with random parameter values, from ranges listed in Table 2. Then, to somehow simulate real measurement, additional white Gaussian noise was added to each signal. Signal to noise ratio was equal to 0 dB for every signal, which means that the power of added noise was equal to the power of the original signal. Next, for each noised signal, the parameters were obtained 10 times using the proposed method.
This process is graphically presented in Figure 4. For clarity, only one signal generation was visualized.
Each reference signal was generated with 1 kHz sampling frequency using a mid-point numerical method of solving differential equations and lasted for 20 s.

Real, Measured Signal
The proposed method was additionally tested on a real signal. To measure it, a 14channel Emotiv EPOC+ NeuroHeadset electroencephalograph was used. Emotiv EPOC+ NeuroHeadset is designed for scalable and contextual human brain research. Arrangement of specific electrodes utilize 10-20 standard. The signal was measured between the O2 electrode and the Common Mode Sense (CMS) electrode [32], which was chosen because alpha waves have biggest amplitude in occipital lobe, especially in the nondominant hemisphere [33]. Signal was sampled with 128 Hz frequency. Obtaining electrode resistance not exceeding 10 k Ohm at the start of the measurement is an essential condition for obtaining a good quality EEG. It depends, to a large extent, on proper preparation of the scalp, which before the electrodes are attached, should be carefully degreased and the superficial calloused layer of the epidermis removed.
The participant was a 25-year male without any diagnosed neurological disorders. The measurement was taken while the subject was well-rested and without drug usage in prior 6 months.
After putting on the electroencephalograph, the participant relaxed, closed their eyes, and after some time, opened their eyes again. To avoid rise of frequency and amplitude of the signal, immediately after closing eyes persisting for a few seconds [34], only part of the signal was used for further analysis (seconds 15-25 in Figure 5a,b). This part, with corresponding discrete Fourier transform (DFT), is shown in Figure 5c,d.
The measurement was taken while the subject was well-rested and without drug usage in prior 6 months.
After putting on the electroencephalograph, the participant relaxed, closed their eyes, and after some time, opened their eyes again. To avoid rise of frequency and amplitude of the signal, immediately after closing eyes persisting for a few seconds [34], only part of the signal was used for further analysis (seconds 15-25 in Figure 5a,b). This part, with corresponding discrete Fourier transform (DFT), is shown in Figure 5c,d. To minimize the impact of artefacts, including the imperfect adhesion of electrodes to the head skin, the acquired signal was first filtered in the range of 2-20 Hz, so range greater than the typical alpha wave frequencies (8)(9)(10)(11)(12)). They are characterized by an amplitude of 3 to about 50 μV on average [36]. The activity of these waves increases when the subject lay with eyes closed before falling asleep. It is also present when a person is at rest, in a relaxed or meditative state, etc. Alpha waves disappear during increased mental effort [37].
During the conducted research, the digital bandpass filter with steepness parameter of 0.85 and stopband attenuation of 60 dB was used.

Results
In Table 3, the results of the repeatability and the accuracy studies are shown.  To minimize the impact of artefacts, including the imperfect adhesion of electrodes to the head skin, the acquired signal was first filtered in the range of 2-20 Hz, so range greater than the typical alpha wave frequencies (8)(9)(10)(11)(12)). They are characterized by an amplitude of 3 to about 50 µV on average [36]. The activity of these waves increases when the subject lay with eyes closed before falling asleep. It is also present when a person is at rest, in a relaxed or meditative state, etc. Alpha waves disappear during increased mental effort [37].
During the conducted research, the digital bandpass filter with steepness parameter of 0.85 and stopband attenuation of 60 dB was used.

Results
In Table 3, the results of the repeatability and the accuracy studies are shown. Parameters obtained for the real signal are shown in Table 4. Search ranges and values used in original work were also included.
The convergence curve of the genetic algorithm is presented in Figure 6. It can be seen that the initial random population had the best individual about 2.4 worse than the final population (2.02 vs. 0.85). Parameters obtained for the real signal are shown in Table 4. Search ranges and values used in original work were also included. The convergence curve of the genetic algorithm is presented in Figure 6. It can be seen that the initial random population had the best individual about 2.4 worse than the final population (2.02 vs. 0.85). Figure 7 shows 5-s fragments of the original signal (7a) and the modeled signal (7c) with corresponding spectra (plots 7b and 7d, respectively). Even though the signals are slightly different in the time plane, their frequency planes are almost identical. It is to be expected, as the function minimized by the algorithm was defined in the spectral plane.  Figure 7 shows 5-s fragments of the original signal (7a) and the modeled signal (7c) with corresponding spectra (plots 7b and 7d, respectively). Even though the signals are slightly different in the time plane, their frequency planes are almost identical. It is to be expected, as the function minimized by the algorithm was defined in the spectral plane. After evaluating defined function for 100,000 times, the mean time needed for one evaluation was determined to be about 0.013 s. Multiplying this value by the number of individuals in each population (256) and the number of populations (150), the maximal time for the whole genetic algorithm was calculated to be about 9 min (Intel Core i7-3740QM CPU with four 2.7 GHz cores, 16 GB RAM). Time is the main disadvantage of the proposed method. It was not focus of this study, however, to obtain results in an online manner. After evaluating defined function for 100,000 times, the mean time needed for one evaluation was determined to be about 0.013 s. Multiplying this value by the number of individuals in each population (256) and the number of populations (150), the maximal time for the whole genetic algorithm was calculated to be about 9 min (Intel Core i7-3740QM CPU with four 2.7 GHz cores, 16 GB RAM). Time is the main disadvantage of the proposed method. It was not focus of this study, however, to obtain results in an online manner.

Discussion and Conclusions
In this paper, a genetic algorithm was used to obtain one column Jansen-Rit model parameters. A spectral similarity function was minimized to generate most optimal parameters for a specific individual. This method proved to be repeatable and accurate, but quite computationally expensive. However, given the convergence curve in Figure 6, the best value was obtained as a best score already in about the 30th population. That means that algorithm could potentially give results faster. On the other hand, convergence of stochastic algorithms tends to depend strongly on initial population. In future works, it would be beneficial to conduct repeatability studies on the proposed method with limited maximal population. Results of conducted repeatability studies with the proposed maximal population suggest that the algorithm does not get stuck in local minima.
Results of repeatability and accuracy can point out which parameters have the greatest contribution to the value of the minimized function, and, therefore, in future studies, which parameters could be expected to constitute the most informative signal features. Hence: A and B (the maximal post-synaptic potential amplitudes) have rather high accuracy compared to the other parameters. ICC value of A, however, is comparatively low. Nonetheless, by itself, the intraclass correlation of 0.897 can be interpreted as good [31]. C, so the numbers of synapses between neuron populations, proved to be the most determinative of the minimized function, wielding ICC at 0.978 and the accuracy of 0.853. It is not surprising, as this parameter is known to have biggest influence on the shape of the modeled signal [10]. Parameters describing the activation function gave mixed results. On one hand, v0 and r turned out to be quite descriptive. On the other hand, ICC and accuracy values of the e0 parameter suggest that it is not as informative. It may be due to the quite narrow sought range and should be studied in future research. Parameters characterizing the noise signal delivered to the model also proved to be very descriptive, wielding highest accuracy and very high ICC. Parameter p was extensively studied in [10], providing deeper insight into its influence on the whole model.
Values of the parameters obtained for the real signal significantly differ from the original works and are close to defined ranges' extremes. That could suggest that those ranges could be further widened in future works.
Results obtained in this paper show that the genetic algorithm works well in a J-R model parameter obtaining task. The defined function can be effectively minimized. Due to the stochastic character of this type of algorithm, the main disadvantage of the proposed method is the computational cost-it is certainly not a fast way to obtain model parameters. However, there is already an online method of parameter obtainment [23][24][25]. State observers proposed in those works predict model parameters dynamically, i.e., in a specific timestamp. The method proposed in this paper enables to obtain parameters for arbitrarily long time periods. Furthermore, this research was conducted on a personal computer. Utilizing some sort for computation cluster would significantly reduce the computation time.
Utilization of J-R model parameters as a classifier input should be checked in future studies. The parameters should be obtained for an already classified set of signals and used to train some classification algorithm. Institutional Review Board Statement: Ethical review and approval were waived for this study, due to not performing any invasive measurements or treatments to a participant, which was also the author of the paper.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.

Conflicts of Interest:
The authors declare no conflict of interest.