It is known that chemical reactions are responsible for acquiring information, transmitting it, and decision-making in living organisms. However, the number of man-written algorithms that can be efficiently executed using a chemical medium Adleman, (1994); Kuhnert (1986), Kuhnert et al. (1989); Steinbock et al. (1995); Agladze et al. (1997); Vazquez-Otero et al. (2012) is quite limited. We believe that the difficulties in practical applications of chemical computers are mainly caused by the fact that an efficient strategy for signal coding using a chemical medium has not been developed yet. The most popular approach to chemical computing imitates information processing with semiconductor devices. In semiconductor devices, information is binary coded in different values of electric potentials.

The success of semiconductor technology came as the consequence of the highly efficient realization of semiconductor binary logic gates. Such gates are characterized by a long time of error-free operation. The semiconductor gates can be assembled together inside an integrated circuit, producing more complicated information processing devices. The bottom-up design strategy Feynman et al. (2000) perfectly matches the technology; complex devices are made as a concatenation of simpler ones. On the contrary, chemical logic gates Toth and Showalter (1995); Steinbock et al. (1996); Sielewiesiuk and Gorecki (2001); Adamatzky et al. (2002) are not small nor fast. The time of their reliable operation is measured in hours, not in years. It seems very hard to make a chemical medium where millions of gates are combined together and work as planned for a long time.

Here we describe and discuss an example of chemical computation based on non-binary information coding. The presented results are based on numerical simulations of the time evolution of the considered computing medium. Such an approach has been motivated by similarities between nerve signals and propagating pulses in a spatially distributed medium in which Belousov-Zhabotinsky (BZ) reaction proceeds Gorecki (2015). The BZ-reaction is catalytic oxidation of an organic substrate in an acidic environment Belousov (1959); Zhabotinsky (1964); Field and Burger (1985); Epstein and Pojman (1994). Two stages of BZ reaction can be visually identified. One of these stages is the fast oxidation of the catalyst. The other is a slow reduction by an organic substrate. The color of the solution of BZ-medium reflects concentrations of catalyst in the oxidized and reduced forms. Therefore, the nonlinear behavior of the medium as oscillations between reduced and oxidized states, propagation of the region characterized by a high concentration of oxidized catalyst, or appearance of spatio-temporal patterns can be easily observed. In a spatially distributed medium, where BZ-reaction proceeds, a local excitation corresponding to the high concentration of HBrO ₂ can propagate in space in the form of the concentration pulse. This type of behavior resembles the propagation of nerve impulses in living organisms. As a result, the BZ-reaction has attracted attention as a medium for experiments with neuron-like chemical computing Adamatzky et al. (2005); Gorecka and Gorecki (2006). Within the most popular approach to chemical computing with BZ-medium, it is assumed that information is transmitted by propagating pulses of the oxidized form of catalyst. For a binary coding, the presence of a pulse represents the logic TRUE state, and the state with a low concentration of the catalyst in the oxidized form is the logic FALSE state Field and Burger (1985); Epstein and Pojman (1994).

If the ruthenium complex ( R u ( b p y ) 3 2 + ) is used as the reaction catalyst, then BZ-reaction becomes photosensitive Kadar et al. (1997). Oscillations can be inhibited by light. For the same initial concentrations of reagents, the medium can oscillate at dark and show an excitable behavior at a low light intensity. And it converges to a steady-state when it is strongly illuminated.

In a medium with photosensitive BZ-reaction excitable channels in which signals can propagate can be formed by specific illumination of a spatially distributed medium. Using a suitable geometry of excitable and non-excitable channels, one can force an appropriate type of interactions between excitations and, for example, make a signal diode Agladze et al. (1996), a memory cell, or logic gates Adamatzky et al. (2005); Yoshikawa et al. (2009). However, in typical applications, such gates are big (with an area of about 1 cm²), and a single operation takes more than 10 s Epstein and Pojman (1994). Therefore, the bottom-up approach from gates to complex information processing tasks does not look promising if the binary is used with BZ-medium.

Recent studies have demonstrated that an oscillating BZ-reaction can be efficiently applied for information processing Smelov and Vanag (2018); Vanag VK. (2019); Proskurkin et al. (2020); Egbert, (2019); Duenas-Diez and Perez-Mercader (2019); Duenas-Diez and Perez-Mercader (2020). For example, it has been shown that a network of interacting chemical oscillators can be trained to perform classification tasks with a reasonable accuracy Gruenert et al. (2015); Gizynski and Gorecki (2017a). To illustrate the problem, let us consider a database D _A composed of N records: D A = p n 1 , p n 2 , … , p n k , q n , n = 1 , N . (1)

The records have a form of (k + 1) tuples ( p n 1 , p n 2 , … , p n k , q n ) , where the first k elements are predictors and the last element (q _n ) is the discrete record type. A classifier of D _A is supposed to return the correct data type if the predictor values are used as the input.

We can easily define the predictor values and the corresponding record type for many life-inspired classification problems, but usually we do not know how to relate both quantities. For example, the needs of medical diagnostic belong to such class of problems. The input information (the predictor values ( p n 1 , p n 2 , … , p n k ) ) are collected from several medical tests on the patient ^#n. On this ground, we are expected to conclude if the patient is healthy or not (the value of q _n ). Our knowledge of the relationship between input (results of medical tests) and output (patient’s health condition) is based on previously accumulated examples. For such problems, the top-down design strategy Gizynski et al. (2016) of a chemical computing medium seems to be more beneficial than the bottom-up one.

Let us assume we have selected a classification problem to be solved. In the following text, it is the determination if a patient has schizophrenia or not. To apply the top-down strategy, we should decide about the medium that is supposed to perform the classification. Here we assume that a network of interacting chemical oscillators can approximately solve the determination of schizophrenia problem. We do not know which network is the best for this task, so we considered a few simple networks illustrated in Figure 2. We selected these particular networks because our previous studies indicated that for a fixed number of modes, a large number of connections in increases the accuracy. All networks can include nodes (oscillators) of two types Gizynski and Gorecki (2017a); Gizynski et al. (2016); Gizynski and Gorecki (2016). The input nodes belong to the first type. They are used to introduce the values of predictors into the network. If a node is assigned as the input of the ith predictor, then its oscillations are suppressed for the time interval of length related to the value of p _i . There are also so-called normal oscillators that are inhibited for a fixed time that is not related to the predictor value. These normal oscillators moderate interactions between oscillators in the medium and optimize the network to solve a specific problem. The time intervals during which their activity is suppressed do not depend on the input. These time intervals define the program executed by the network. For the analysis presented below, we assume that the output information is coded in the number of oscillation cycles observed on a given node. The choice of the output oscillator follows directly from the network optimization. The complete definition of a computing network includes the number of oscillators in the network, their types, locations, the information about the time intervals they are active, and the information about interactions between oscillators.

Obviously, a network with randomly selected parameters has a small chance to work as a good classifier. We have to optimize its parameters (i.e., to teach a network) to perform the selected function. Teaching means that we need a teacher, and in our optimization, it is a specific database T _A that contains diagnostic results and information if a patient is ill or healthy EEG (n. d.). In the following, we do not change the number of oscillators in a network nor modify the geometry of interactions between them.

The application of the top-down strategy to the considered networks means that the parameters such as locations of the input and normal oscillators, inhibition times for the normal oscillators, the method for inputting the values of predictors, or the parameters of reactions responsible for interactions between oscillators are the subjects of optimization. The optimization is supposed to achieve the best match with a representative (training) dataset of cases T _A . We have found Gruenert et al. (2015); Gizynski and Gorecki (2017a); Gizynski et al. (2016); Gizynski and Gorecki (2016) that evolutionary optimization oriented on obtaining the best classifier for a representative training dataset of the problem can lead to a computing network that performs the anticipated task with reasonable accuracy.

In previous papers on chemical database classifiers Gruenert et al. (2015); Gizynski and Gorecki (2017a); Gizynski et al. (2016); Gizynski and Gorecki (2016) an oversimplified event-based-model reflecting the basic features of the oscillator time evolution and of interactions between oscillators coupled by mutual activations was used. The event-based-model divides an oscillation cycle into three phases: excitation, refractory and responsive phase. It also assumes a sharp difference between these phases. An oscillator in the refractory phase is not susceptible to stimulations by interacting oscillations. However, the event-based-model allows for the excitation of an oscillator in the responsive phase that is in contact with an excited oscillator. In this paper, we consider a more realistic model. We represent the time evolution of an individual oscillator using the two-variable Oregonator model I. R. Epstein and Pojman (1994); Field and Noyes (1974) of the photosensitive Belousov-Zhabotinsky (BZ) reaction.

If we neglect interactions with the other oscillators of the network, then equations describing the time evolution of jth oscillator are: d u j d t = 1 ε ( u j − u j 2 − f v j + ϕ j t u j − q u j + q ) (2) d v j d t = u j − v j (3)

Where the variables u _j and v _j represent concentrations of an activator (U _j ) and an inhibitor (V _j ) for proceeding reactions. The parameter ɛ sets up the ratio of time scales for variables u and v, q is a scaling constant, and f is the stoichiometric coefficient. We used the same values of model parameters for all oscillators in the network in our simulations: ɛ = 0.2, q = 0.0002, and f = 1.1. The parameters of the Oregonator model were fixed and did not undergo optimization.

If we assume that the time evolution of oscillators is described by a model of photosensitive BZ-reaction, then oscillators can be individually controlled by illumination, and we include this feature into the considered model as the time-dependent function ϕ _j (t) in Eq. 2. The time-dependent function ϕ _j (t) that describes the influence of illumination on an oscillator is proportional to the light intensity. We considered ϕ _j (t) in the form: ϕ j t = 0.1 ⋅ 1.001 + tanh − 10 ⋅ t − t i l l u m j (4)

In this definition t _illum (j) > 0 defines illumination of the jth oscillator. In the time interval [0, t _illum (j) − δ] (δ = 0.1) the value of ϕ _j (t) is high ( ∼ 0.2 ) . The Oregonator model with parameters given above predicts a stable steady state corresponding to u _j = 0.0002 and v _j = 0.0002. For long times (t > t _illum (j) + δ) the value of ϕ _j (t) approaches 0.0001 what corresponds to oscillations with the period of approximately 10.8 time units. The value of δ describing the speed of transition between the steady state and oscillations can be reduced by increasing the multiplier under tanh () function. In order to extract the answer of a classifier we consider the number of activator maxima within the time interval Z = [0, t _max]. It can be noticed that Eq. (4) has a physical meaning for any value of t _illum (j). If t _illum (i) < 0 then ϕ _i (t) is small and the oscillator ^# i is active during the whole observation interval. When t _illum (k) > t _max than ϕ _k (t) is large and the oscillator ^# k is blocked within Z and does not oscillate. Figure 1 shows the time evolution of activator u ₁(t) predicted by Eqs 2 and 3 with t _max = 99.3 and t _illum (1) = 28.4 (cf. Table 3). As seen for such parameters, oscillations restart just after t _illum (1), and the system produces seven activator maxima within the observation time interval.

In the investigated networks the values of t _illum (j) for normal oscillators were fixed. On the other hand, if an oscillator was considered as the input of p ⁱ (i = 1, 2) then t _illum (j) was an affine function of p ⁱ (cf. Eq. 7). Such function is defined by two parameters t _start and t _end and it has the form: t i l l u m j = t s t a r t + t e n d − t s t a r t ⋅ p i (5)

We assumed that t _start and t _end are the same for all predictors in the considered schizophrenia records (p ¹ and p ²).

It was demonstrated that even a small network composed of 16 or fewer oscillators, with the time evolution of mutual excitations described by the event-based-model could be used to diagnose if a cancer cell is malignant or benign Gizynski and Gorecki (2017a). In this report, we concentrate on designing a network of oscillators that can determine whether the patient has schizophrenia. Schizophrenia is one of the most common forms of psychotic behavior. The patients experience hallucinations, delusion, or disorganized speech. However, schizophrenia is difficult to detect Siebenhuhner et al. (2013). It is believed that the analysis of EEG signals recording brain activity can help to verify if a patient is ill or healthy Ben et al. (2007). The relevant EEG signals were recorded from electrodes placed in different parts of the scalp (see Figure 3). We postulate that a network of interacting chemical oscillators in the form presented in Figure 2 can process the information extracted from the EEG signals and help diagnose schizophrenia. Preliminary results obtained using the network illustrated in Figure 2D were reported in the extended abstract of the ASPAI 2020 Conference Bose and Gorecki (2020). For the database T _A , available on the web EEG (n. d.) containing signals recorded on N = 84 patients, out of which N _h = 39 were healthy and the other had symptoms of schizophrenia (N _s = 45) the optimized classifier returned 82% correct answers. The extended study on the schizophrenia classifier in the network’s pentagon geometry (Figure 2D) was published in the International Journal of Unconventional Computing (IJUC) Bose and Gorecki (2021). In this IJUC paper, we also studied if the classification accuracy can be improved by dividing the whole recorded signal into 3 shorter (20 s long) parts and processing these shorter signals separately. Separated networks with Figure 2D geometry were optimized for each time subinterval within 500 generations of evolution. Next, the majority procedure was applied to obtain the final classification results. Such a method increased the classification accuracy of records in T _A to 90%. In the current report, we investigate if the classification accuracy can be improved by the network geometry. The results of schizophrenia diagnosis using classifiers with geometries illustrated in Figures 2A–C are new and have not been previously reported. Moreover, the presented results for the pentagon geometry differ from those published because the genetic optimization was performed for 260 more steps than in Bose and Gorecki (2021) and a new maximum of fitness was achieved during these additional steps. Therefore the pentagon-shaped classifier has a different structure than previously reported.

The manuscript is organized as follows. Transformation of EEG Signals Into the Input Data describes how the input data are extracted from the EEG signals. In Numerical Model of Information Processing Network, we present a numerical model for the simulation of network time evolution. Network Optimization and Results gives details of network optimization. The conclusions and suggestions for the future development of the networks for schizophrenia diagnosis are presented in the following Conclusion and Discussion.

The considered networks were small, and there was no room for too many input variables because each predictor requires its input oscillator. We used the signals recorded from F7 and F8 channels marked red in Figure 3 as the inputs for schizophrenia detecting networks described below. Such choice is motivated by previous studies indicating that the signals obtained from the frontal lobe of the brain reveal the difference in the brain activity between a schizophrenic patient and a healthy subject Ben et al. (2007).

The EEG signals were recorded with a sampling rate of 128 Hz for 1 min. The medical EEG test of each patient produced 16 data files corresponding to signals recorded on different electrodes. Each data file contained K = 7,680 values of recorded potential (in μV). The time between consecutive potential values is Δt = 7.812 5 ms. Let V ^l (n, k) denote the potential recorded for nth patient, on the lth electrode and at the time t _k = k ⋅Δt. To reduce the size of input data, we averaged the recorded signals. Therefore, each signal was trimmed to a single number. We assumed that time average signals provide us with a sufficient amount of information to diagnose schizophrenia. The averaged potentials were defined as: x n l = ∑ k = 0 K V l n , k (6)

Next, the time averaged potentials recorded for the whole set of patients were normalized. We introduced: μ l = 1 N ∑ n = 1 N x n l and σ l = 1 N − 1 ∑ i = 1 N x n l − μ l 2

The values of predictors p n 1 and p n 2 for the patient n were defined as: p n 1 = x n F 7 − μ F 7 σ F 7 , p n 2 = x n F 8 − μ F 8 σ F 8 (7) In the above equations μ ^F7 = 7.724 μV, μ ^F8 = 2.46 μV, σ ^F7 = 20.3 μV and σ ^F8 = 15.10 μV.

As a result, the problem of schizophrenia diagnosis is reduced to the best classification of the training dataset: T S = { ( p n 1 , p n 2 , q n ) , n = 1 , N } where the record type q _n = 0 for a schizophrenic patient and q _n = 1 for a healthy subject. The distribution of records in the T _s database in the (p ₁, p ₂) coordinates is illustrated in Figure 4. Blue and red crosses correspond to schizophrenic and healthy cases, respectively. It can be seen that the points corresponding to different cases are not separated, which makes their classification difficult.

In our study, we also considered predictors generated using combinations of signals recorded on other electrodes, but for those, the classification accuracies on corresponding datasets were lower.

The time evolution of reactions proceeding in a single oscillator was described by the two-variable Oregonator model (Eqs 2 and 3). We assumed that interactions between the oscillators ^#k and ^#j appeared as the results of reactions involving the activators U _k and U _j of these oscillators: U j + B j → U k + C k (8) U k + B k → U j + C j (9) with identical reaction rate constant k _B .

We also assumed that the activator of each reaction could spontaneously decay in the binary reaction: U j + D j → p r o d u c t s (10) with the rate constant k _D . In the equations above symbols B, C and D denote other molecules involved in these reactions. Reactions (8–10) involving species B, C and D are formally introduced to explain chemical communication between nodes and to justify the mathematical description of it.

Therefore, the changes in concentrations of U _k and U _j as the result of reactions (8 and 9) are: d u j d t = − k B b j u j (11) d u k d t = − k B b k u k (12) and the changes in concentration of U _j as the result of reaction (10) is: d u j d t = − k D d j u j (13) In Eqs 11–13 b _j , b _k and d _j denote concentrations of B _j , B _k and D _j , respectively. We assume that these concentrations were high with respect to concentrations of activators involved and the same for all oscillators. Therefore, the concentrations of B, C, and D were regarded as constant during the network evolution, and there is no need to include them in the model of network evolution. Let us introduce symbols α and β defined as: α = k _D d _j and β = k _B b _j . Keeping in mind that values of α and β can be modified by concentrations of B _j and D _j , we can treat them as free parameters that can be easily adjusted. Therefore, the values of α and β can be included in the optimization procedure. The same mathematical description of interactions between nodes applies for controlled exchange of reaction mixtures between nodes and outflow of activator in a system with immobilized catalyst.

Within our model the time evolution of the network is described by the following set of kinetic equations: d u j d t = 1 ε u j − u j 2 − f v j + ϕ j t u j − q u j + q − ( α + β ∑ i = 1 , m s j , i ) u j + β ( ∑ i = 1 , m s j , i u i ) (14) d v j d t = u j − v j (15) The last two terms in Eq. 14 represent the coupling in between ith and jth oscillators and the activator decay. The symbols s _j,i are defined as:

s _j,i = 0 if j = i or if j ≠ i and oscillators ^# j and ^# i do not interact,

s _j,i = 1 if j ≠ i and oscillators ^# j and ^# i do interact.

The set of Eqs 14 and 15 describes the network evolution after all parameters characterizing the medium including t _illum for all oscillators are known.

A classifier is supposed to produce an answer within a finite time. However we do not know it. Therefore, the time t _max that defines the interval of time for which the network evolution is observed [Z = (0, t _max)] is one of the optimized parameters of a classifier. We postulate that information about patient health is extracted from the number of activator maxima recorded on a selected oscillator of the network, during the time interval Z. In order to find which oscillator should be used as the output one we calculate the mutual information I (G; O _j ) Cover and Thomas (2006) between the discrete random variable G of record types in the training dataset T _S (G = {q _n , n = 1, N}) and the discrete random variable O _j of the number of activator u _j maxima o _j (n) observed on the jth oscillator in the network when the predictors of nth database record are used as the network input (O _j = {o _j (n), n = 1, N}). The mutual information I (G; O _j ) can be calculated as: I G ; O j = H G + H O j − H G , O j (16) where H () is the Shannon information entropy Shannon (1948) and the discrete random variable (G, O _j ) = {(q _n , o _j (n)), n = 1, N}. The oscillator ^# i for which the mutual information between G and O _i is maximal is used as the network output. The mutual information calculated for the output oscillator was considered as the measure of network fitness: F i t n e s s = m a x j I G ; O j (17) It can be expected that in the majority of cases the optimization based on the mutual information leads to a classifier with the highest accuracy Gorecki (2020).

One of the most significant challenges of civilization is how to use Artificial Intelligence (AI) for various life-inspired problems AI (2022). AI techniques can be beneficial for medical applications where the knowledge is accumulated as information on previously cured cases. To diagnose a new patient, one should search for similarities with the previous ones. Here we presented an application of AI methods for designing a system that can help to diagnose schizophrenia. We assumed that schizophrenia could be detected by a chemical oscillator network that analyses EEG signals recorded from electrodes located on a patient scalp.

We considered a few information processing networks characterized by different numbers of nodes (c.f. Figure 2). We think that networks of interacting chemical oscillators represent more realistic models of biological neural computing than typical artificial neural networks with arbitrarily selected activity rules Kay (2003). Their time evolution is described by realistic kinetic equations that model specific nonlinear chemical reactions. Here we used the Oregonator model for the Belousov-Zhabotinsky reaction. We believe that other, more realistic models Vanag and Epstein (2009) of chemical evolution lead to qualitatively similar results. It seems that oscillator networks require a smaller number of nodes than standard neural networks to achieve the same accuracy. The results presented in Gorecki and Bose (2020) indicate that just three oscillators can solve a geometrical problem of how a point from a unit square is located with respect to a disk placed at the square center with 95% accuracy. The fact that just a few oscillators can perform a complex information processing function, confirmed by the results for schizophrenia diagnostics reported in this paper, opens the door for experimental realization of chemical instant machines with systems of interacting oscillators reported in the literature Juan Manuel et al. (2020); Vanag and Yasuk (2018); Vanag. (2019a); Proskurkin et al. (2019); Gizynski and Gorecki (2017b).

We expected that the schizophrenia diagnosis accuracy increases with the network size. However, the networks formed of 3, 5, and 6 oscillators gave 82% of correct answers for cases included in the training dataset. On the other hand, this number is much higher than given by standard classification methods included as options in the Clasiffy procedure of the Mathematica program Mathematica (2021). If we apply this procedure to the training dataset T _S , the highest classification accuracy (76.2%) is obtained for the GradientBoostedTrees method, whereas the NeuralNetwork option leads to 73.8% accuracy. The accuracy of schizophrenia diagnosis using chemical oscillator networks increased to 86.9% if three networks process a case and the majority rule is used to select the final answer. The fact that larger networks did not produce better results than the small ones could be related to inefficient optimization for a large number of parameters that were taken into account. The problem can be overcome by a larger population of classifiers and a larger number of optimization steps. However, both methods increase the numerical complexity of the optimization. It is also worth mentioning that the recent study on the application of machine learning methods for schizophrenia detection from textual input Wawer et al. (2022) reported ∼ 80% accuracy on a sample of 94 people (47 ill and 47 healthy).

Although the presented results are encouraging, datasets with a larger number of patient data are important for further studies. A dataset with a large number of patient records can be separated into a training dataset of a few thousand cases and much larger testing ones Gorecki and Bose (2020) that is independent of the training one. The separation of records between the training and the testing datasets can be done in many different ways. By selecting different training datasets, one can verify the stability of the schizophrenia diagnostic classifier with respect to different training. The observation that classifiers with similar parameters are obtained for different training datasets confirms that the diagnostics are unbiased by selecting a training dataset.

The presented classification method is based on many assumptions. All of them can be lifted in search of the best network for the schizophrenia diagnosis.

The optimization of interactions in the medium can be directly included in the optimization program. The presented results were obtained assuming that interactions between oscillators were fixed, as shown in Figure 2. The information about interactions was included in equations describing the time evolution of oscillators as the binary parameters s _i,j in Eq. 14. The values of these parameters can be included in classifier optimization. The network model includes the activatory coupling between oscillators. It means that an excited oscillator can speed up the excitations of the other oscillators that interact with it. Such coupling is observed, for example, in droplets containing reagents of BZ reaction. Alternatively, one can consider a medium with inhibitory coupling where excitation of one oscillator slows down the activity of those oscillators that interact with it Vanag and Yasuk (2018); Vanag. (2019b); Proskurkin et al. (2019). Allowing for different types of coupling within a single network can help to identify the best medium for a given computing task.

We assumed that the output could be related to the number of activator maxima observed at a specific oscillator. However, one can consider alternative methods of extracting information from network evolution Borst and Theunissen (1999); Zhang et al. (2014). For example, the output can be related to a pair of numbers of activator excitations recorded on two selected oscillators Gorecki (2020).

It is anticipated that the accuracy of diagnosis should improve if the information on signals recorded on more than two electrodes is included in the input. The presented optimization algorithm can be easily modified to do this if one includes input oscillators of any important signal into the network. If additional signals do not increase the Fitness, then networks with the inputs of irrelevant signals will vanish from the population. Moreover, to improve the accuracy of a large oscillator network, one should consider different decay rates α at different oscillators and different activator transfer rates β for individual couplings. Of course, it results in a significant increase in the number of parameters undergoing optimization.

Future studies should reveal if the generalizations of classifiers as described above can significantly increase their accuracy of schizophrenia diagnosis if compared with the classifiers presented in this paper. After a successful classifier optimization is completed, its application does not require significant computing power. There are just two steps of the algorithm: 1) normalization of patient data with parameters (μ, σ) obtained for the training dataset and 2) numerical solution of differential Eqs 14 and 15 and activator maxima counting. A modern laptop needs a few seconds to execute these tasks. The whole procedure can be incorporated into EEG equipment software or distributed as a laptop or smartphone application.