Design of MRI structured spiking neural networks and learning algorithms for personalized modelling, analysis, and prediction of EEG signals

This paper proposes a novel method and algorithms for the design of MRI structured personalized 3D spiking neural network models (MRI-SNN) for a better analysis, modeling, and prediction of EEG signals. It proposes a novel gradient-descent learning algorithm integrated with a spike-time-dependent-plasticity algorithm. The models capture informative personal patterns of interaction between EEG channels, contrary to single EEG signal modeling methods or to spike-based approaches which do not use personal MRI data to pre-structure a model. The proposed models can not only learn and model accurately measured EEG data, but they can also predict signals at 3D model locations that correspond to non-monitored brain areas, e.g. other EEG channels, from where data has not been collected. This is the first study in this respect. As an illustration of the method, personalized MRI-SNN models are created and tested on EEG data from two subjects. The models result in better prediction accuracy and a better understanding of the personalized EEG signals than traditional methods due to the MRI and EEG information integration. The models are interpretable and facilitate a better understanding of related brain processes. This approach can be applied for personalized modeling, analysis, and prediction of EEG signals across brain studies such as the study and prediction of epilepsy, peri-perceptual brain activities, brain-computer interfaces, and others.

www.nature.com/scientificreports/ classification module. In the first step, the bio-signal is encoded as a spike train using spike encoding algorithms and then, the generated spike sequences are entered into the SNNr for unsupervised learning using the spike-time dependent plasticity rule (STDP) 15 . The outputs of the SNNr neurons are connected to a SNN-based classifier/ regressor trained in a supervised mode 11,21 . The connection weights in the SNNr are updated according to an unsupervised STDP learning algorithm, which makes the network more biologically plausible, but the quality of its performance depends strongly on the encoding stage and the initial connection weights in the SNNr that are set following a small-world connectivity model 15,21,25 . The paper addresses some new challenges in this field, namely: • Instead of using general brain templates, we propose to use a personal MRI data for pre-structuring a personal SNNr model, MRI-SNNr, that could lead to a better analysis of personalized EEG data and better accuracy of EEG signal prediction. • Developing new learning algorithms for EEG signals in a MRI-SNNr, based on real value EEG data, avoiding encoding the data into spike sequences. • Using the models from above to predict brain signals in areas not measured and not used for training the models. • Developing new methods for the analysis and interpretation of trained MRI-SNNr for a better understanding of the data and the personal brain processes.
The proposed here 3D MRI-SNNr architecture uses the Izhikevich neuronal model 16 allowing for continuous EEG data to be processed in the SNNr and to avoid a preliminary spike encoding process. Using the Izhikevich neuron requires solving two differential equations governing the membrane processes and membrane potential of neurons in a model. These equations are: where: v is the membrane potential; I denotes the input current to the neuron as a continuous value; u represents recovery current; a, b, c and d are constants, affecting the neuronal spike behavior. The voltage is in millivolts, and the time is in milliseconds.
The article is organized as follows: The second section introduces the proposed method and MRI-SNNr architecture for personalised modelling along with the learning algorithms. "Results and discussions" presents simulation results of the performance of two MRI-SNNr models when compared with the existing NeuCube model and with univariate EEG signal prediction method. "Discussions and Conclusion" suggests possible applications of the proposed method and directions for further research.

The proposed method and MRI-SNNr architecture for personalised modelling and learning algorithms
The method and the overall architecture. The proposed MRI-SNNr architecture consists of a SNNr and Output module (Fig. 1). The SNNr is structured according to MRI personal data, consisting of interconnected observed (input/output) and hidden neurons. The boundary element method (BEM) 26 using FieldTrip tool 27 and MATLAB software to localize brain activity are used to structure the 3D-SNNr using a personal MRI image data. The 3D locations of the extracted points from the Calculated-Volume-mesh grid are assigned to the spatial positions of spiking neurons in a 3D SNNr. The observed neurons are positioned proportional to the spatial locations of the EEG channels providing input information 15 . Figure 2 presents the algorithm for the spatial localization of the spiking neurons in the 3D SNNr following spatial information from both MRI and EEG data.
After structuring the MRI-SNNr, EEG data are used to train it. The connection weights between the spiking neurons are first initialized using small world connectivity rule 15,21,25 as it is expected a higher impact of neighboring neurons so as to reduce the deviation of the normalized output EEG signal from the real EEG signal, described in the next section. In this paper, neural synaptic connection weights of the observed neurons are updated based on gradient descent learning rule.
The regulator/classifier module has the task of normalizing and filtering (0 ~ 60 or 70 Hz) 28,29 the output data from the observed neurons to produce continuous value output signals (Fig. 1).
In this framework, it is presumed that there is a bi-directional link between neurons in a neighborhood, and each neuron output is connected to the same neuron as a feedback. N i denotes the neighborhood of i th neuron as follows: where T is the maximum spatial distance of two neurons that are chosen depending on the problem. Here T = 70 mm, and D ij is the distance between neuron i and j. In the case of i = j, feedback from the output of the neuron to its input occurs. www.nature.com/scientificreports/ Similar to 9 initial bi-directional connection weights w ij = w ji = 0 between neighboring neurons i and j in the defined above SNNr are set before training of the SNNr with EEG data (see also Eq. 2): The connection weights of the observed (input/output) neurons are modified with the use of the presented below new GDR learning algorithm. The connection weights of the hidden neurons are modified during learning of EEG data using the STDP rule and after learning, only the strongest connection between every two neurons is kept as suggested in 9 .
(3) w ij = wji = 0 j ∈ N i no connection j / ∈ N i Figure 1. A schematic diagram of the proposed MRI-SNNr architecture for EEG signal modelling and prediction using personalized MRI data for structuring the SNNr ("hidden" neurons, in yellow,) and EEG input data for training the SNNr (input/output neurons, also called observed neurons, in red). The proposed learning algorithms for the MRI-SNNr. We propose that the connection weights of the observed neurons are modified using the gradient descent rule (GDR) as described further below and the connection weights of the hidden neurons are modified with the use of the spike-time dependent plasticity (STDP) learning rule. The key motive of the GDR and STDP combination is to reduce the dependency of SNN on spike encoding algorithms, as GDR uses real value data, and to preserve biological plausibility still using STDP for spiking information processing in the hidden neurons. For completeness, both variants of continuous and spike inputs to the models are implemented and their performance is compared as in some studies only spike information is available.
The hidden neurons are updated according to the STDP learning rule 15,21 : where w ij is the synaptic weight change in N i and Δt is the time difference between the spike time of the ith neuron and the postsynaptic neuron j. A + and Aare the maximum values of synaptic modification in the two directions and τ + and τrepresent the time constant of permissible changes in synaptic weights. The method for the modification of the observed-neuron-synaptic weights is presented below. With regard to Eq. (1) and the relationship between membrane potential and excitation input current, this paper proposes a new formula for the calculation of the input current I i to the ith observed neuron of the Izhikevich type: where w ij is the synaptic weight between neurons i and j and dθ i dt is the leakage current expression of the ith neuron. In fact, the suggested Eq. (5) can mitigate the calculation cost of the error backpropagation process ahead, while the effect of the neural activity of adjacent neurons is considered. In the other word, this paper has proposed Eq. (5) as the input current adaptation law with the aim of creating of linear algebraic summation relationship among membrane potentials of neurons in a neighborhood so as to reduce the cost of gradient descent calculation, while preserving the biological concept of the neural network by means of utilizing Izhikevich-type neurons.
The relationships for the observed neurons are rewritten in matrices form as follows: where the parameters are described as below: and if j / ∈ N i then w ij = 0 , where N is the total number of reservoir neurons in the SNNr and N ob is the set of observed neurons.
A more general form of Eq. (1) is considered as the following: Then, considering G(V ) = g(V )dt , F(v) = Idt = tanh(WV) + θ , on case of zero initial conditions, the neuronal membrane potential is calculated as: The cost function J (t) at time t is defined as: where e i (t) is the normalized ith neuron output error at time t, denoted by v n i = v i (t) , from the output of the EEG channel corresponding to the location of that neuron, represented by v d i (t).
Substitute Eq. (11) into Eq. (10), the following equation is obtained: Accordingly, weight modifications for the observed spiking neurons are determined based on the GDR: Thus, the rule of updating the observed neuronal synaptic weights concerning the effect of neighboring neurons is as follows: where T s is the simulation time step. The leakage current expression in the update law is calculated as: Concerning Eqs. (4) and (14), combination of STDP and GDR, as a reinforcement learning strategy, is described as: where u s (−�t) is a step function, defined as: Two realisations of the proposed MRI-SNNr architecture. The proposed MRI-SNNr architecture and learning algorithms are utilized here in two realizations: (a) continuous value EEG signals are used as inputs (MRI-cSNNr); (b) spike sequences are used as inputs as in many studies only spike sequence data is available (MRI-sSNNr) (Fig. 3).
In both realizations from Fig. 3, neurons in the MRI-SNNr are located using the same personalized MRI data. Parameters of the neurons and STDP learning are initialized according to Table 1. In the MRI-sSNNr, every spike sequence X is converted into a continuous value signal using a baseline B(t) with regard to the method from 9 and with regard to the raw input signal S(t), before the data is used for training the SNNr applying the same learning algorithm as in 2.2.
With regard to Eq. (1), and with the aim of real neural activity simulation corresponding to the input EEG signals, it is suggested that the general form of observed neuron's potential be considered according to the desired potential and the actual EEG value of the corresponding channel (Eq. 19).
where X is obtained according to the following equations: where V EEG is the real input EEG signal. After calculating

Results and discussions
To compare the performance of the proposed in "The proposed method and MRI-SNNr architecture for personalised modelling and learning algorithms" two MRI-SNNr models with other methods for EEG signal modelling and prediction, two cases of study data have been used. In this regard, "Experimental results of the proposed

STDP connection weights updating parameters
Parameter Value Parameter Value www.nature.com/scientificreports/ personalized MRI-SNNr models compared with other EEG predictive modelling methods on the first case study EEG data", "Would adding more neurons in the MRI-SNNr result in an increased accuracy of EEG signal prediction?", "Can a MRI-SNNr model be used to predict signals at other locations of the SNNr corresponding to non-monitored brain areas?", "Analysis of the experimental results" and "Neuronal activity analysis for a better understanding of brain activities" present analysis of experiments on the EEG data of the first case study, and "Experimental results on a second case study MRI and EEG data" presents experiments on the second case study data. The first case study data were collected from the Epilepsy Long-term EEG Monitoring Center of the University Hospital (Imam Khomeini), approved by the Research Ethics Committee of the Neuroscience Research Institute of the Tehran University of Medical Sciences. It is carried out as part of a clinical study with the relevant guidelines and regulations, and the patient has given informed consent to this. The second dataset was collected from the freely accessible database http:// eeg. pl/ epi, and the experiment has been conducted during a routine procedure at the Warsaw Memorial Child Hospital. The patients with refractory epilepsy were selected for a pre-surgical study with the aim of removing the epileptogenic zone.
Patients characteristics. Patient  Experimental results of the proposed personalized MRI-SNNr models compared with other EEG predictive modelling methods on the first case study EEG data. Here, we compare the performance of the proposed two MRI-SNNr models from Fig. 3 with two other methods. The first one is a traditional univariate (a single EEG signal prediction) method, also called "random walk" that predicts the next time moment signal to be equal to the previous. The second model for comparison is NeuCube 15 , which uses spike encoded EEG signals for learning and the SNNr is structured using a general Talairach brain template rather than using a personalized MRI data. The NeuCube model has 1471 neurons in the SNNr, located at brain areas according to the used template 15 . The NeuCube software from (http:// www. kedri. aut. ac. nz/ neucu be) is used.
The experimented MRI-SNNr has 300 Izhikevich neurons consisting of 15 observed (EEG channels) and 285 hidden neurons. Modeling properties of the network are set such that 20% of neurons are inhibitory and the rest are excitatory. In this paper, we considered that the time unit of Ts = 8 ms is appropriate to reconstruct the lower -than-70-Hz-frequency data, which includes all spike-based frequency bands. For a plausible comparison of the outputs of the models, both the input EEG signals and the network outputs are filtered simultaneously using a MATLAB-designed-Minimum-Ordered-Equripple-FIR filter (MOEF) with a 0-70 Hz pass frequency and 80 dB attenuation high frequency, containing the main frequency bands and eliminating drift (all filter parameters are listed in the Appendix) . The model parameter values used for the experiments are presented in Table 1.
To evaluate the performance of the proposed MRI-SNNr models and to compare it with other models, simulation is performed with the same EEG data. 70% of the data is used for training and 30% for testing.
To evaluate the performance of the proposed MRI-SNNr models and to compare it with other models, simulation is performed with the same EEG data consisting of 500 samples measured every 8 ms. 70% of the data is used for training (350 samples, 2.8 s) and the future 30% for testing (150 samples, 1.2 s, but in Figs. 4 and 6 only the first 0.2 s are visualized). During training, input data are entered one by one and the observed neurons are trained to predict their next time output (in 8 ms) by using the output error and the GDR algorithm. During training, hidden neurons update their synaptic weights by using STDP law. During testing, there is no update of the connection weights of observed neurons. During the test procedure, no data samples are entered and the model predicts 1.2 s of the signal ahead. In this stage, only STDP is utilized in order to update the weights in the hidden neurons, unsupervised. In this regard, the observed neuron's output is measured and compared with the actual values to calculate the test error and all the indicated parameters in Table 2 has been calculated based on all 500 samples error.
Prediction results for the two MRI-SNNr models are presented in Fig. 4a for each of the 15 EEG input channels, with 0-70 Hz output MOEF filter applied on both the outputs and on the input data for comparison of results. In order to work out the effect of the filter characteristics, Kaiser filter with 0-70 Hz frequency band and 80 dB attenuation at higher frequencies, has been imposed in Fig. 4b as the post-processor of raw EEG signals and also post-processor of outputs of the suggested network. A glance at Fig. 4 reveals that the shape of filtered outputs and inputs is vastly relevant to the sort of utilized filter and the other characteristics of the output filter such as frequency passband. The proposed continuous value-based MRI-cSNNr outperforms the spike-based MRI-sSNNr in terms of convergence, speed and prediction error in both a and b. With regard to the lower prediction error while utilizing the MOEF 0-70 Hz, subsequent simulations are performed with this filter. In a following section we will explain how such models can predict outputs at other neurons which correspond to unmonitored/ not measured area of the brain (Fig. 4c).  Table 2 for both structures. The estimated MSE for the 300-neuron SNNr models is not substantially different from the 100-neuron network MSE, while the estimated MSE for the 400-neuron SNNr models has dropped significantly. While increasing the number of neurons is not expected to have a linear effect, more neurons allow for more areas of the brain to be investigated. The determination of the appropriate values for the neuron parameters can influence the prediction accuracy, although the broad analysis of this issue is not the aim of this paper. The next simulations therefore proceed with 300 neurons and the same initialization.

b). (c)
Predicted signals at 10 neurons that correspond to unmonitored signals not used for training, while the models were trained on the measured EEG data from the 15 channels. It is seen that the prediction accuracy within the MRI-cSNNr structure is higher. www.nature.com/scientificreports/

Can a MRI-SNNr model be used to predict signals at other locations of the SNNr corresponding to non-monitored brain areas?
As it is well accepted in the literature for focal epilepsy localization for example, the more EEG channel information is used, the higher the accuracy of focal localization is achieved 30,31 . On the other hand, there can be limitations in recording channel information and some hospitals can monitor only a limited number of channels, say around 32, simultaneously. Would a 3D MRI structured SNN model provide a possibility of predicting signals at other locations of the model in addition to the EEG channel locations (observed neurons)? If this is achieved, it may become possible to detect changes in other parts of the brain where no EEG channels are located. And this will be due to the use of MRI personal information to structure the SNNr. In this section the performance of the two MRI-SNNr models from Fig. 3 to predict signals at SNNr locations that correspond to non-recoded EEG channel is evaluated. For this purpose, 15 channels of EEG signals are considered as inputs and 25 EEG channel locations in the 3D SNNr model are examined as outputs, including the additional 10 locations of hidden neurons that correspond to presumed non-observed 10 EEG channels. Figure 4c illustrates 10 EEG channels prediction performance to validate the estimation of non-recorded EEGs. The 10 continuous value signals, not used for training but evaluated for prediction, correspond to EEG channel positions O2, O1, P4, P3, C4, C3, F3, F4, Fp1, and Fp2 not used in the training data. MSE of the EEG signal prediction across all 10 new sites as well as 15 inputs' prediction is shown in Table 2.
Proper convergence in the prediction of un-observed sources indicates the potential of this method in predicting unmonitored brain areas which constitutes the first study in this respect. Table 2 presents some statistical analysis, entailing maximum Standard Deviation (SD) of prediction error of all predicted channels and the total mean square error (MSE) 31 of EEG prediction of several methods when compared with the proposed MRI-SNNr methods for 0-60 Hz and 0-70 Hz output module filtering. The table shows that the proposed method for EEG signal modelling and prediction based on the proposed MRI-SNNr approach has resulted in several orders of magnitude better prediction accuracy. The following equation describes the MSE calculation:

Analysis of the experimental results.
where E is the normalized-one-step-ahead-prediction error vector, outlined within Eq. (10).
In order to create a suitable comparison between NeuCube and the novel proposed methods, the F3 channel is also considered as the unmonitored signal output of NeuCube. The experiment indicated that NeuCube has the capacity of predicting unmonitored signals as much as the suggested structure, but its' prediction error is higher.
Another set of experiments was conducted with 100 and 400 spiking neurons instead of 300 as reported in Fig. 4. Results of all experiments above are presented and validated in Table 2.
It can be seen from Table 2 that selection of the output filter parameters, such as type of the filter, stop band, pass band and attenuation, influence the accuracy of the models. Lowest error was registered for 400 spiking neurons in the SNNr when compared with 100 and 300. P-Value analysis, as another validation criterion, also indicates that both real data and predicted ones, are statistically analogous.
Neuronal activity analysis for a better understanding of brain activities. Connection weights analysis of the MRI-cSNNr model could be used to better understand functions of the brain. In our model, the activation degree of each neuron is defined as suggest in 9 by: www.nature.com/scientificreports/ where w ij and w ji are the weights of bidirectional connections between neurons i and j which included negative and positive values. Figure 5 shows the degree of activation of each of the 3D neurons in the SNNr alongside the neuronal potential distribution and its firing pattern in 11, 78 and 135 epochs of incremental training of the MRI-SNNr models. epochs.
Here the EEG signals were reported from a patient undergoing hyperventilation (HV). The pattern of degree of activation converges after 100 epochs to a specific pattern, which can be further explored in future research when compared with a pattern produced by standard stimuli. Having studied some researches in the field of HV effects on epileptic EEG neural activity rhythms 4,5 and owing to what can be inferred here, the degree of activation of the brain's central and parietal lobes under hyperventilation rises compared to other regions and the prediction of the signals in all measured regions have been accurately achieved as presented in the previous sub-sections.
Experimental results on a second case study MRI and EEG data. In order to further validate the proposed method, EEG signal prediction test is developed for the second case study data in a similar way as for the first case study, with 15 monitored inputs, and 32 outputs, including 17 non-monitored channels (Fig. 6). The non-monitored channels are P4, T6, O1, O2, S1 ~ S11, A1, A2. Figure 6 reveals that the suggested strategy is successfully able to predict EEG signals in both status, namely, 70% training and 30% test of monitored data or 100% testing of non-monitored EEG channels as a part of total data. Figure 7 illustrates the similarity between the power spectrum of the input EEG FP1 signal (signal 1), after being filtered for comparison with the signal from the output module, also filtered. We can see that a similar power spectrum of the signal is produced by the model to the one of the input EEG signal which is also an indication for an accurate modelling result.
As it is shown in the above figures, although the signals in Fig. 4 and Fig. 6 may visually look similar to a sine wave, the frequency analysis shows that the recorded signals contain more widespread frequency bands rather than a mere sinusoidal wave. The wave signals for all channels in all experiments in this paper are visualized after applying the same filtering. Some of the wave forms may look similar to a sine form when perceived visually, but the spectrum analysis shows a complex frequency spectrums rather than a single frequency.

Discussions and conclusion
In this paper, a new method is introduced for the creation of personalized SNN models, MRI-SNNr, for the analysis, prediction and understanding of EEG signals, where MRI data of a person is used to structure spatially a 3D SNN model. The 3D SNNr model consists of observed (input/output) neurons and hidden neurons. A new machine-learning algorithm is proposed for the MRI-SNNr to learn the EEG data using gradient descent learning rule for the observed Izhikevich neurons and STDP learning rule for the hidden neurons. Both continuous value inputs and spike inputs are considered in two separate versions of the proposed MRI-SNN architecture. www.nature.com/scientificreports/ The presented preliminary experimental results on case study personal EEG data manifest a better performance of the proposed method applied on continuous value EEG signals rather than on spike encoded signals, both resulting in several orders of magnitude better modelling accuracy than a traditional statistical method and the NeuCube model 15 , which is structured according to a general brain template rather than a personal MRI information.
The learned neural connections in the MRI-SNNr, when using input data to the observed neurons to train a model, makes it possible to predict signals in other neurons of the model that correspond to non-measured (not observed) brain areas. This method can be interpreted for a better understanding of larger scale brain activity across applications. The strong performance of the MRI-SNNr in the prediction and reconstruction of signals www.nature.com/scientificreports/ of the SNNr corresponding to unmonitored EEG sources (data not used for training the models) suggests that this approach can be potentially used also for localizing critical brain regions. This is the world first research that demonstrates that a MRI pre-structured personalised SNN model trained on data from a number of EEG channels, can be used to predict the activity of other areas of the brain corresponding to EEG channels that are not used for training the model. The proposed method was illustrated on predictive learning and analysis of EEG data of the epileptic patients. Dynamic changes in EEG spatial-temporal patterns in different brain areas were discovered, which could be explained by transitions between epileptiform, preictal, and ictal events 28 . For the used case study of personal MRI and EEG data, it was found through analysis of the MRI-SNNr model that the right temporal region has the highest potential, which is consistent with the reported location of focal epilepsy in that area for this subject. The method was also illustrated on a second case study MRI and EEG data, and the high potential of hidden states of the brain functions would stem from the ability of 17 other non-monitored source prediction. The proposed method has the potential to be used for further study and prediction of epilepsy events 32,33 .This requires further investigation with the use of data from other individuals. The proposed MRI-SNNr architecture and learning algorithms are of a generic type and could be potentially used across wide range of brain studies, such as peri-perceptual studies 34 , brain-computer interfaces 22 and other 21 .
The accuracy of modeling depends on many parameters, related but not limited to, how MRI-SNNr is structured to map a personal MRI; what are the optimal parameters of learning in the SNNr; what are the optimal output functions and filter parameters. More research is planned to examine and optimize these parameters for a better personalized predictive modelling of EEG data. Representing spatio-temporal patterns of activities in the MRI-SNNr during learning in the form of spatio-temporal symbolic rules is also a challenging task that is also considered for future research by the team 13,21 .
Received: 12 July 2020; Accepted: 9 April 2021 Figure 7. Power spectra of the input EEG Fp1 signal (signal 1) and the output signal from the corresponding neuron in the SNNr model, both filtered in the same way for a proper comparison of the output.