EEG channels selection for stroke patients rehabilitation using equilibrium optimizer

2.1 Inspiration of EO

The EO draws its inspiration from the concept of mass balance applied to a control volume. This approach involves using the equation of mass balance to account for the concentration of nonreactive constituents within the control volume. The mass balance equation ensures that the inflow, outflow, and generation of mass within the control volume are taken into consideration and maintained. Several parameters need to be taken into account when formulating the mass balance equation, including the rate of change of mass over time, the amount of mass entering the control volume, the amount of mass generated within the control volume, and the amount of mass leaving the control volume. The formulation of the mass balance equation is as follows:

where V is the control volume, C presents the concentration amount in V , Q is the volumetric flow rate in and out of V , G is the mass generation rate in V , and C eq denotes the equilibrium state concentration.

To obtain the concentration in the control volume, equation (1) can be rearranged in the following manner:

subject to:

By incorporating integration, equation (2) can be expressed as follows:

This results in:

Subject to:

The concentration in control volume V , given a known turnover rate, can be estimated using equation (4), where C 0 represents the initial concentration at time t 0 . Furthermore, equation (4) can also be utilized to calculate the average turnover rate when the generation rate is known.

The aforementioned equations play a crucial role in controlling the search behavior of the EO. In the EO population, each particle corresponds to a solution, and the concentrations represent their positions (analogous to swarm algorithms). As described earlier in equation (4), the update mechanism for a particle’s position and concentration involves three distinct terms: the equilibrium concentration, the concentration difference between the particle and the equilibrium state, and the generation rate. These terms collectively govern the adjustment of the particle’s position and concentration during the optimization process.

2.2 Procedure of EO

The EO is a metaheuristic algorithm used for solving optimization problems. Here’s a summary of the procedure based on the article you provided:

2.3 Fitness function evaluation

During this step, all particles within the population undergo evaluation using the fitness function specific to the optimization problem at hand. Subsequently, the particles are sorted based on their corresponding fitness values. This sorting process allows for the identification of the equilibrium candidates among the particles, as those with better fitness values are considered more favorable solutions.

2.4 Equilibrium pool and candidates

Once the EO search process begins, equilibrium candidates are identified to establish the search pattern for the particles. Through extensive testing across different problem types, it has been determined that the best four equilibrium candidates exhibit superior performance during the optimization process. In addition, a fifth candidate is introduced, which represents the arithmetic mean of the four best candidates. This inclusion aims to strike a balance between exploration and exploitation capabilities within the EO. These five candidates are utilized to construct the equilibrium pool vector, as depicted below:

During each iteration of the EO, all particles update their concentrations by randomly selecting candidates with equal probabilities. For instance, in one iteration, a particle may update its concentrations based on C → eq ( ave ) , while in the next iteration, the concentrations could be updated using C → eq ( 1 ) . This process continues, with particles updating their concentrations using different candidates in a randomized manner, until the stop criterion is met. The stop criterion signifies the end of the optimization process.

2.5 Exponential term ( F )

The exponential term ( F ) is employed to update the concentrations in the EO. F plays a crucial role in enhancing and achieving an optimal balance between exploration and exploitation during the search process. The calculation of F can be performed using the following formula:

where λ → is a random vector within a range between 0 and 1 and t represents an iteration function that can be computed as follows:

where Iter and Iter Max denote the current and the maximum number of iterations, respectively, and a 2 represents a static value. Equation (9) is utilized in the EO to enhance its convergence while simultaneously maintaining a balance between exploration and exploitation.

where a 1 is a static value. When a 1 is set to a higher value, it increases the exploration capability of the EO, resulting in a lower exploitation capability. Conversely, when a 2 is set to a higher value, it enhances the exploitation capability while reducing the exploration capability. The term sign ( r → − 0.5 ) is employed to influence the directions of exploration and exploitation. Here, r → represents a random vector ranging from 0 to 1. The sign function is applied to the difference between r → and 0.5, determining whether it is positive or negative. This sign determines the direction in which the exploration or exploitation will occur. By utilizing equation (9) for equation (7), equation (10) is extracted.

2.6 Generation rate ( G )

The generation rate is a crucial parameter in the proposed algorithm, as it significantly impacts the quality of solutions by enhancing the exploitation phase. In the context of engineering applications, various strategies can be employed to represent the generation rate as a function of time. One such strategy involves utilizing a first-order exponential decay process to represent the generation rate [22]. The following equation defines this representation:

where G ( t ) represents the generation rate at time t , G 0 denotes the initial generation rate, and λ represents the decay rate parameter. By incorporating this exponential decay process, the generation rate can be appropriately modeled over time. In the EO, the equations employed to reformulate the generation rate are formulated as follows:

Subject to:

where r 1 and r 2 represent random numbers generated between 0 and 1. The generation rate control parameter (GCP) vector, obtained from equation (14), is used to compute the GCP variable. The GCP variable contributes to the updating process, influencing the generation rate control parameter.

A considerable number of particles in the EO utilize the generation concept to update their states. The selection of these particles is determined by a term called generation probability (GP), which is calculated using equations (13) and (14). The decision-making process based on equation (14) depends on the particle’s level. If the GCP is zero, the G is assigned zero, and all the decision variables of these particles are updated without involving G . To strike a balance between exploration and exploitation, the GP is set to 0.5. This value ensures a good trade-off between exploring new solutions and exploiting promising ones. The updating rule of EO is applied as follows:

As indicated in equation (15), the updating process for the concentrations of the particles in the EO involves three terms:

1. The first term corresponds to the equilibrium concentration, which is randomly selected from the elite equilibrium pool. This pool comprises the top five solutions or equilibria obtained during the optimization process.

2. The second term serves as an exploratory operator that helps guide the particles toward the optimal point in the search space. It contributes to the exploration of the solution space.

3. The third term acts as an exploitation operator that enhances the quality of solutions. It is primarily controlled by the generation rate (as described in equation (12)). The generation rate influences the trade-off between exploration and exploitation.

The global and local search behavior is mainly governed by several parameters, including the concentrations in particles and equilibrium candidates, the turnover rate ( λ ), and the sign relationship between the second and third terms. When the signs of these terms are the same, it promotes global search by increasing the diversity. On the other hand, when the signs are opposite, it facilitates local search by reducing the diversity. By adjusting these parameters, the EO can strike a balance between global and local search, leading to effective exploration and exploitation in the optimization process. Figure 1 shows the movements of the equilibrium candidates.

Figure 1

Equilibrium candidates’ collaboration in updating a particles’ concentration [21]. Source: Created by the authors.

Although the second term in equation (15) aims to target solutions from search space regions that are distant from the equilibrium candidates, and the third term extensively explores the search space region that includes the equilibrium candidates, and it is important to note that this behavior is not always guaranteed. Small changes made to the denominator of the third term can result in increased variations, which can contribute to the exploration of certain dimensions. The second term in equation (15), represented by C 1 − C eq , and the third term, represented by C eq − λ ⋅ C 1 , are influenced by the generation rate terms described in equations (13) and (14). The variation in these terms is governed by the value of the generation rate parameter λ , which can vary depending on the dimension.

Interestingly, high variation is more pronounced in dimensions where λ has smaller values. This mechanism is analogous to the mutation operator commonly observed in evolutionary-inspired algorithms. It allows for increased exploration in dimensions that may require more attention or have not been adequately explored, promoting the search for potentially better solutions in the optimization process.

2.7 Particle’s memory saving

To enhance the performance of the EO and enable particles to track their trajectory in the search space, memory-saving procedures have been incorporated. In each iteration, the fitness value of a particle is compared to its fitness value in the previous iteration. If the current fitness value is better, it replaces the previous value, thus preserving the particle’s best position so far. This mechanism reinforces the exploitation capability of the EO, as particles are able to retain their best-known solution. However, it is important to note that this mechanism can also introduce a potential challenge of stagnation in local minima if the algorithm does not perform sufficient global search. In such cases, the EO may become trapped in local optima and struggle to escape to potentially better solutions in other areas of the search space.

2.8 Exploration ability of EO

The exploitation task in the EO is achieved through the utilization of various parameters and mechanisms, including:

• a 1 is a parameter in the EO algorithm that plays a crucial role in controlling the exploration capability. It determines the distance between a new position and the equilibrium candidate. A higher value of a 1 increases the exploration ability of the EO algorithm.

• sign ( r − 0.5 ) is employed to control the trajectory of exploration. The variable r is generated randomly with a uniform distribution between 0 and 1.

• GP parameter utilizes the generation rate as a probability for contributing to the concentration updating process. The GP value determines the extent to which the generation rate term influences the exploration performance during optimization. When GP is set to one, it indicates that the generation rate term has a high contribution to the exploration performance. Conversely, when GP is set to zero, it implies that the generation rate term will always be contributing to the optimization process. This increases the likelihood of the algorithm getting stuck in local optima, as exploitation is prioritized over exploration.

• The equilibrium pool consists of five particles. The selection of these five particles is based on a slightly random criterion, which has been determined through initial experiments.

2.9 Exploitation ability of EO

The exploitation task is obtained utilizing the following parameters and mechanisms:

• a 2 is responsible for the exploitation task, specifically in determining the extent of exploitation. It governs the amount of deep search conducted around the best solution found so far.

• sign ( r − 0.5 ) controls the direction of exploitation and drives the local search in a specific direction. The variable r is randomly generated from a uniform distribution between 0 and 1.

• Memory saving is a memory parameter that keeps track of the best particles found so far during the optimization process.

• Equilibrium pool at the end of iterations. The equilibrium candidates are becoming more closer to each other.

The EO pseudocode in Algorithm 1 shows a comprehensive overview of the EO steps.

3.1 Phase I: EEG signal acquisition

The EEG data of poststroke patients with hemiparesis of the upper extremities was studied. This study involved three poststroke patients treated with the recoveriX system (g.tec medical engineering GmbH), with a mean age of 22 years ( SD = 4.582 ). Each participant got BCI-based MI training for three months, with two training sessions each week (for a total of 25 training sessions). The study team conducted and analyzed two assessments (pre- and posttraining). The pretraining evaluation was scheduled 30–35 days prior to the intervention, whereas the post-training evaluation was conducted a few days after the intervention (see Figure 2). The Ethikkommission des Landes Oberosterreich in Austria ( # D − 42 − 17 ) authorized this study protocol, and each patient signed an informed consent form prior to the preassessment. Finally, this dataset is captured with a sample rate of 256 Hz. According to system indications for a MI mental task, the patients were instructed to visualize dorsal wrist movement. Before and after EEG neurorehabilitation, the patients participated in 25 MI-based BCI sessions with follow-up evaluation visits to measure the functional changes. Each session consisted of 240 MI repetitions with both hands and broken down into three 80-trial runs. Each session lasted approximately 1 h, including the time required for setup and cleanup. The MI-based BCI tasks were illustrated with randomized inter-trial intervals in a pseudo-random order.

Utilized were EEG caps with sixteen active electrodes from g.Nautilus PRO, g.tec medical engineering GmbH, Austria, based on the international 10–20 system.

3.2 Phase II: Preprocessing

Each channel of the recorded EEG dataset initially used two standard filters. First, a bandpass filter (BPF) with frequencies between (8 and 30 Hz) was used to confine the band of the recorded EEG dataset, [23], and second, a butter-worth ( B W ) notch filter at (50 Hz) was used to reduce the noise from power line interference. ( I C s ) s ( t ) = [ s 1 ( t ) , … , s n ( t ) ] utilizing the F a s t I C A algorithm proposed by [24]. The set s ( t ) of n unknown components that were linearly mixed by within matrix A and x ( t ) is the set of n observations where x ( t ) = [ x 1 ( t ) , … , x n ( t ) ] [25–27], represents the EEGs and are related to s ( t ) , t is the time, the ICA equation is calculated as follows:

then, three metrics were used to evaluate the artifactual components ICs: Kurtosis (Kurt), skewness (Skw), and sample entropy (SampEn). If these parameters surpassed the ± 1.2 threshold for each IC, the IC was marked as critical and denoised using WT. The practical value of the threshold is determined through trial and error and previous research [28–30]. The threshold value of ± 1.2 is not a drawback of AICA-WT technique as the artifactual ICs will not be rejected but will be denoised using the WT technique. As a result, WT denoised the marked ICs before returning the enhanced components to the original EEG dataset [31,32].

3.3 Phase III: EEG features extraction

Three different feature extraction techniques are used in this work. These techniques include time domain features, entropy features, and frequency domain features as extracted in [7].

3.4 Phase IV: EEG channel selection using EO optimizer

This phase is the main contribution of this work, and it includes several steps to achieve the optimal subset of the EEG channel, which can provide the highest accuracy rate. The following steps represent how can we adapt the EO optimizer for the EEG channels selection problem.

• Step 1: EO parameters initialization for EEG channel selection. In this step, the EO parameters are initialized with initial values. These parameters are the number of microbats (solutions) ( N ), maximum number of iterations ( Max Itr ), max ( C max ) and min ( C min ) frequency range, and the bandwidth range ( ε ) in the range [ − 1 , 1 ].

Also in this step, the EO population memory for the EEG channel selection problem should be Initialized. The EO solutions are randomly generated in this step using equation (17), considering the EO algorithm and particular optimization problem parameters and constraints.

where ∀ i = 1 , 2 , … , d d refers to the total number of the EEG channels, ∀ j = 1 , 2 , … , N , and r i is a random number between 0 and 1. The produced solutions together generate the population, as shown in equation (19) and are stored in the EO memory (EOM) in the ascending order on the basis of their fitness values. The best solution with the fittest values is assigned to EEGsol Gbest and stored in the first position of EEGEOM.

• Step 2: EO population intensification. Now, all microbats EEGsol will change their position considering by a frequency f generated randomly, as shown in equations (20) and (21). Accordingly, the EO positions will be updated using equation (22).

Notable that the microbats are updating their locations to be closer to the global best ( EEGsol Gbest ). Thus, the new position of the microbats intensifies the position with a direction to EEGsol Gbest .

• Step 3: EO population diversification.

In this step, the microbats positions are updated based on random parameters in attempting to find better global solutions. A new solution is selected from the BM based on a selection method and assigned to EEGsol best . Subsequently, the current solution is updated using the EEGsol best , as follows:

where A ′ is the loudness average. As Steps 2 and 3 are two contradictory steps, EO algorithm will choose one of them for execution at a time. This selection is based on the following equation:

• Step 4: EEGBM update. The position of the microbats in the EEGBM will be updated in this step if the new location is fitter than the old one and R ≤ A j . Moreover, EEGsol Gbest will be updated if the new BM contains a solution with a better fitness value. Subsequently, r j and A j values will be updated in accordance with equations (25) and (26).

where itr denotes the current iteration number.

• Step 5: Stop criterion. Steps 2, 3, and 4 will be repeated until the algorithm reach the stop criterion.

3.5 Phase V: Classification and performance measures

In this work, the k-nearest neighbor classifier is used. The performance of the proposed method and the optimal solution of the EO optimizer for the EEG channel selection problem was evaluated using the values of average classification accuracy. The overall average accuracy is computed using the following equation:

where x refers to the number of correctly classified input data points and N refers to the total number of attributes.