Developing A Seizure Prediction Algorithm for A Non-Invasive Neuromodulator

In this study, a novel real-time seizure prediction algorithm is introduced to predict epileptic seizures. The proposed algorithm is expected to be applicable in a noninvasive neuromodulator. As a model of the epileptogenic zone, a small-world network of Huber-Braun neurons was built up. To assess the effects of noninvasive stimulation techniques, such as transcranial magnetic stimulation, this network was modified and the magneto-motive forces and the electromagnetically induced currents were further applied on the network. Comprehensive investigations of the electroencephalograms of epilepsy patients have suggested that the seizures are generated by some chaotic mechanisms. Hence, chaos and bifurcation theory were applied and the induced current was considered as the bifurcation parameter. The bifurcation diagram of the 'inter-spike' intervals of the mean voltage of the small world network was obtained. The precise time at which the bifurcation took place was subsequently considered as the time of the seizure onset. Comparisons of the bifurcation diagrams obtained from the patients’ electroencephalographs showed that the proposed network model can reasonably represent the actual neuronal networks of the epileptogenic zone. A dataset of the electroencephalographs of epilepsy patients and normal volunteers from an epilepsy center in Germany was used to validate the prediction algorithm. The simulation results show that the proposed algorithm has a significant capability to predict the precise occurrence of seizures and the achieved sensitivity, accuracy, and specificity of this approach were remarkably higher than those reported in previous studies.

 Abstract-In this study, a novel real-time seizure prediction algorithm is introduced to predict epileptic seizures. The proposed algorithm is expected to be applicable in a noninvasive neuromodulator. As a model of the epileptogenic zone, a small-world network of Huber-Braun neurons was built up. To assess the effects of noninvasive stimulation techniques, such as transcranial magnetic stimulation, this network was modified and the magneto-motive forces and the electromagnetically induced currents were further applied on the network. Comprehensive investigations of the electroencephalograms of epilepsy patients have suggested that the seizures are generated by some chaotic mechanisms. Hence, chaos and bifurcation theory were applied and the induced current was considered as the bifurcation parameter. The bifurcation diagram of the 'inter-spike' intervals of the mean voltage of the small world network was obtained. The precise time at which the bifurcation took place was subsequently considered as the time of the seizure onset. Comparisons of the bifurcation diagrams obtained from the patients' electroencephalographs showed that the proposed network model can reasonably represent the actual neuronal networks of the epileptogenic zone. A dataset of the electroencephalographs of epilepsy patients and normal volunteers from an epilepsy center in Germany was used to validate the prediction algorithm. The simulation results show that the proposed algorithm has a significant capability to predict the precise occurrence of seizures and the achieved sensitivity, accuracy, and specificity of this approach were remarkably higher than those reported in previous studies.

I. INTRODUCTION
Epilepsy is one of the chronic neurological disorders that affects 1% of the world population. This disease consists of a set of abnormal excessive synchronized discharges that arise in some cortical neurons [1]. In 25% of epilepsy cases, seizures are resistant to various conventional treatments including anti-convulsive drugs and epilepsy surgery. This type of epilepsy is known in the literature as refractory epilepsy [2]. Alternative seizure control protocols have been introduced for controlling seizures in these cases. These protocols usually have two steps. The first is in the form of detection or a prediction algorithm (PA) and the other is to design a neuro-stimulation technique. Various researches have been undertaken on the subject matter under investigation, most of which focused on the seizure onset paradigms aimed at developing a PA.
The process in the prediction algorithms involves extracting the features that are related to seizure onset from the online EEG records. It can also be achieved from the responses of the neural models of the epileptogenic zone (EZ). Dynamical studies on the electroencephalograms (EEGs) of epilepsy patients have revealed that the epileptic seizures are generated by some deterministic nonlinear mechanisms. There is strong evidence to suggest that such mechanisms are chaotic and even the occurrence of a minute variation in their initial conditions, can bring about considerable changes in their outputs [3]. A thorough investigation of the EEG records using the modern nonlinear theory (chaos and bifurcation) and extracting the features that are related to chaotic systems, including the correlation dimension and the Largest Lyapunov exponent (LLE), have lead researchers to get novel and deeper insights into the mechanisms of initiating and evolving the epileptic seizures [4] [5]. Moreover, bifurcation theory had widespread use for pattern recognition for several neurological disorders [6] [7][8] [9]. This theory deals with changes in qualitative behavior or the topological structure of a given family. One example of this is the alteration of one parameter in a set of differential equations. The parameter in question is referred here as the bifurcation parameter. The bifurcation theory has been applied for pattern recognition from the EEGs of the pre-ictal phase of epilepsy patients [6]. It has also been applied to achieve the same results from the responses of the in-vitro epileptic neuronal models [7][8] [9].
As pointed out earlier, two major steps involving in controlling the refractory epilepsy are developing an appropriate PA and stimulating the epilepsy focus to mitigate the effects of seizures. The neuro-stimulation therapies are classified into two main classes, namely responsive and non-responsive. Depending on the stimulation technique, each class can be divided into invasive and non-invasive categories. In responsive neurostimulators, the stimulus is triggered after the detection of the seizure onset. However, in non-responsive ones, bursts of stimulus are applied on the EZ, over a predetermined time interval. This renders the seizure prediction unnecessary [10]. The responsive neurostimulation (RNS) [1], nonresponsive Vagus nerve stimulation (VNS) [11] and deep brain stimulation (DBS) [12], can claim to be a few examples of invasive stimulation techniques proposed so far. Their painless mechanism and insignificant side effects are among the major characteristics of the non-invasive neuro-stimulation techniques including transcranial magnetic stimulation (TMS) and repetitive transcranial Developing A Seizure Prediction Algorithm for A Non-Invasive Neuromodulator Mahnaz Asgharpour, Mehdi Sedighi, and Mohammad Reza Jahed Motlagh stimulation (rTMS) that have led to their universal application in the scientific and medical domains. These techniques have therefore been experimentally used by the experts as an effective treatment strategy to manage the seizures in epilepsy patients [3] [13]. The antiepileptic effect of the low-frequency high-intensity rTMS has been recognized [14]. Application of the rTMS on the EZ, while has shown to decrease the frequency of seizures, it has also proved to have a simultaneous ability in improving the psychological status of the epilepsy patients [15]. Despite considerable advances in designing the non-invasive neurostimulators, their sensitivity and accuracy have not been acceptable for clinical applications [3]. The underlying reason being the lack of experimental models to examine the prediction algorithms. Therefore, appropriate modeling of epileptic neurons' behavior, simulating the effects of the non-invasive stimulation techniques on those models, and introducing proper features are vital issues in designing the advanced neuromodulators which should be addressed. Hodgkin-Huxley (HH) model is one of the basic neural models which describes the propagation of action potentials in neurons. Despite the significant effects of the HH-type models in following the biological concept, solving their equations requires great computing power. Hence, they are not appropriate for simulating large-scale neural networks [16]. The Huber-Braun (HB) model of thermally-sensitive hypothalamic neurons, is a simplified HH-type model. It can reproduce a wide range of neural activity, including tonicfiring, bursting and chaotic spiking patterns of the neurons involved in certain cases of neurological disorders such as epilepsy and Parkinson's disease [17] [18]. This model reveals different bursts of impulses as a function of temperature [19] or electromagnetically-induced potentials [8].
In this study, a small-world network (SWN) of the Huber-Braun neurons was presented as a model of the EZ. Owning to the painless and non-invasive nature of magnetic stimulation, the induced effects of this technique on the HB model were simulated. Applying the modern nonlinear theory, some features related to the seizure onset were extracted from the network response, and a PA was proposed. A set of 180 EEG records related to ictal and preictal phases of epilepsy patients as well as normal volunteers from an epilepsy research center in Germany [20], was used for validating the algorithm.
The paper is arranged as follows, the second section deals with the proposed modified small-world-network of the HB neurons. The results and discussions are presented in the third section and finally, the conclusion is presented in the fourth section.

II. METHODS AND MATERIALS
Designing a TMS-based neuromodulator, necessitated simulation of the effects of TMS on a model of the EZ. Given the existence of both local and non-local interconnections in actual neural networks, SWNs was considered to be the most appropriate model for simulating the EZ [18] [21]. Intended to simulate the effects of the electromagnetically induced current, , on an SWN of Huber-Braun neurons, the coil of the stimulator was simulated and the distribution of induced electric field and the activating function was computed. Following this, the cable model of a Huber-Braun nerve fiber was simulated.
The assumption was that a finite myelinated HB nerve fiber is located parallel to the axis of a circular coil, in a homogenous and isotropic volume conductor. The Huber-Braun equations, as well as the methods used for performing the aforementioned simulations, have been presented in the Appendix. This modified HB model was used to simulate the nodes of the network. The chemical synapses were modeled as a means of simulating the connections between the physically far apart neurons.
The membrane potential is the most important variable to investigate different behavior of the classic HB model. The neural messages are transformed via the spikes of the membrane voltage. The inter-spike-intervals (ISIs) that are the distances between two consequent spikes are considered to be the most important variables for calculation [3][22] [23]. The effect of temperature on the HB model behavior has been investigated [19]. In other studies, the Homoclinic bifurcation in the ISIs of the model's membrane voltage has been shown to arise from the variation of its temperature [24] and the injected external current [25]. Moreover, a new bifurcation parameter that is the external electromagnetically-induced voltage, has been introduced in another study [8].
Applying the modern nonlinear theory and considering the induced current as the bifurcation parameter, the bifurcation was detected in the inter-spike-intervals (ISIs) of the mean voltage of the network, (see Appendix). With the early detection of seizure-like activity in the response of an HB neural network in mind, the detection of the bifurcation onset in the ISIs of , at different measures of was assumed as the mane-stone of the prediction algorithm. This is based on the findings of several studies that observed some kinds of bifurcation in the Poincare map of the ISIs of intracranial-EEGs (icEEGs) of epilepsy patients in the pre-ictal phase. For instance, perioddoubling and the Flip bifurcation have been observed in one research on the icEEGs of the patients suffering from temporal lobe epilepsy [22]. Some efforts have however been made to classify the EEGs of the epilepsy patients into ictal and pre-ictal phases with the application of the chaos theory. In one specific study, the LLEs of the ISIs of electroencephalograms were successfully used to classify the different seizure phases with a high degree of accuracy [3]. In this study, the bifurcation was observed in the ISIs of of an SWN that was affected by TMS. Since the bifurcation is a paradigm that has been detected in the preictal phase, the detection of bifurcation in an EEG, could be considered as an indication of the patient experiencing the pre-ictal phase and a seizure is about to emerge. Furthermore, since the bifurcation parameter in this paper is related to TMS, any alteration in that parameter is expected to produce an alteration in the model behavior and control the bifurcation (the seizure onset) in the model.
As a means of validating the results of the proposed prediction algorithm, real EEG datasets are required. The EEG recordings from a database of an epilepsy research center in Bonn Hospital were used for this investigation [20]. These consisted of 120 icEEG recordings related to the ictal and pre-ictal phases of epilepsy patients. There were also 60 scalp-EEGs of 5 normal volunteers. Each recording consists of 23.6 seconds of the EEG data. Owing to the lack of available EEG data under TMS, the simulated effects of TMS was applied on the recordings. It was assumed that a TMS is applied on the EZ, while the EEG is being recorded by a TMS-compatible-EEG-Device. Since TMS-therapy has shown to have anti-convulsive effects that decrease the excitability of the EZ neurons [14][15], the proposed PA was considered an active PA. This meant that the algorithm can use for both seizure suppression and prediction. To put it differently, while the applied nonresponsive magnetic stimulation decreases seizure frequency, the PA predicts seizures that have not been controlled by the nonresponsive stimulation. A stimulus with a different strength or frequency can then be applied to mitigate seizures.
Assuming that , same as the injected electrical currents, has the same value all over the EZ [18], the distribution of the induced voltage, , in a synapse with a 330 ohm-cm resistance was calculated by its impedance relation with the . The EEG signals are not sensitive to the action potentials produced by the cell membrane. It is a generally held belief among scientists that the postsynaptic currents play a major role in producing the EEG signals. After the elimination of noise and artifacts, the EEGrecording-devices amplify these currents and illustrate them in the form of voltage [26]. On the other hand, TMSinduced currents, can depolarize the membrane potential and open the voltage-dependent ion channels to initiate action potentials. However, the synaptic activity subsequent to this phenomenon (the postsynaptic currents) is reflected directly in the EEG. The EEG-recording-device records a linear projection of the distribution of the postsynaptic currents on its measuring channels. EEG signals represent the potential difference of the two scalp electrodes to record an EEG. In the presence of TMS, EEGs from different electrodes are time-locked averaged with respect to the TMS and this averaged EEG signal is the resulted EEG recording [26]. Experimental pieces of evidence suggest that noise and artifact elimination in low-frequency TMS causes a deflection in the EEG baseline [27]. Bearing this in mind, the induced voltages on the EEG-recording-electrodes due to , were amplified and added to the EEG data. It was then assumed that these induced voltages are sampled with the same frequency of the sampled EEG signals. Thereby, applying TMS on EEG signals was simulated.
To analyze the real data, some features such as ISI should be extracted. In this regard, the offline EEG text files were converted into numerical vectors. Intended to compute LLE, the vectors of the sampled EEGs by = 173.61 Hz, were considered as time-series. Using the method developed by Rosenstein [28], the related LLEs were calculated. For each EEG dataset, the bifurcation diagram of ISIs with respect to was presented. In this regard, the estimated values of was added to the real data vectors. The resulted data was divided into distinct 200-sample windows with a duration of 1.15 seconds. Consideration of the prediction speed made it necessary to divide each data vector into as short a window as possible in such a manner to make a prompt detection of the bifurcation possible. The distinct windows should not be too short, because sufficient data might be lost. On the other hand, they should not be too long, as that might affect the prediction speed. By setting the window-length to 200 samples, it was possible to secure feature extraction and bifurcation detection from almost all the EEG recordings. To detect the peaks of EEG spikes, a threshold should be defined so that when a spike passes this threshold, it could be considered as a peak. Regarding the peak-to-peak voltage of each recording, a specific threshold voltage, ℎ was specified to each EEG series. The average time to select ℎ for each recording was approximately 15 minutes.
Intended to design a real-time PA, it was assumed that a 1 Hz TMS is applying to the EZ while EEG is being recorded (it has proven that at the frequencies between 0.9 to 6 Hz, TMS has inhibitory effects on the seizures [14]). Quite similar to [20], it was assumed that the online EEGs are recorded by a 128-channels amplifying system and after passing through a 12-bit analog-to-digital converter, they are recorded by a sampling rate of 173.6 Hz. After investigating each patient's offline EEGs, different patient-specific thresholds were defined as a vector of possible ℎ measures. Each sample was compared to each element of the possible-thresholds vector and the EEG peaks were determined. Then, the ISIs regarding each were computed and the bifurcation diagram of the ISIs versus was extracted. When a bifurcation onset was detected, the algorithm automatically illustrated the related element of the thresholds vector as the patient-specific threshold voltage. After this training phase, the PA could be applied to the EEG samples. Doing this, the real-time computation was performed on each 200-sample vector, then the time difference on the consequent peaks and their related values were registered in other vectors. So, the bifurcation diagram of every distinct window could be drawn and bifurcation onset could be detected.

III. RESULTS AND DISCUSSION
As discussed in the Appendix, the circular coil of a magnetic stimulator with 15 wire-turns and a radius of 1 cm was simulated. It was assumed that the amplitude and frequency of the current passing through the coil being 3800A and 1 Hz, respectively. Then the distribution of the magnetic potential, the induced electrical field, and the activating function were achieved and presented in Appendix. As the next step, a small-world-network of 2000 neurons was constructed. The different behavior of of this network related to the was investigated. The bifurcation diagram of the ISIs of was then extracted (Fig.1). As it has been shown, the SWN, which is a model of EZ under TMS, experiences bifurcation as well as chaotic behavior with respect to . It showed that bifurcation occurs in the pre-ictal phase, so detecting bifurcation in the EEG of an epilepsy patient, could be a cue that the patient is in the preictal phase and a probable seizure is upcoming. To validate the PA, some statistical analysis was performed on the real EEGs and the seizure onset-related features were extracted from the data. Regarding this, the EEGs which were in the form of text-files were converted into numerical vectors. In Fig.2, EEG dataset#1, has been illustrated for different classes of data. As can be seen, the peak-to-peak voltage and therefore the ℎ , are different in the pre-ictal and ictal phases and for the normal volunteer too. Choosing ℎ = 35 by trial and error, ISIs were calculated and illustrated in Fig.3. Conspicuously, the distribution of ISIs around 100 µs is different in each case. Fig.4 illustrates the Histogram of ISIs. For the pre-ictal phase, the histogram has a Poisson-like distribution in comparison with the ictal phase. It has been shown that dynamical reconstructing of ISI in a 2-dimensional state space, could yield efficient information about the dynamics of the main system [23]. So, the dynamical analysis of the recursive maps of ISIs was carried out. The function, , in the recursive equation, ISI +1 = (ISI ), is unknown, but the local behavior of the ISIs' dynamics could be estimated from the Fig.5. The red lines, with the equation, = , have been drawn for better visual analysis.   In the pre-ictal phase, the concentration of the ISI dots along the red line was more than in other states. This can be interpreted to mean that the two consequent ISIs are equal in this phase. Furthermore, there are distinctive differences in the distribution of ISIs in other cases. Moreover, some statistical parameters such as variance, standard deviation, and the peak-to-peak voltage of these EEG signals were computed and presented in table I. The table shows considerable differences between the extracted statistics, such as variance and the peak-to-peak voltage of the preictal and ictal phases. The above computations were performed on other datasets as well. The results revealed that in the pre-ictal phase, these statistics, particularly variance, were dramatically lower than those in the ictal phase. Hence, it seems that a patient-specific algorithm could be defined to classify different ictal phases using these statistics. Applying TMS on the EEG-recordings, a comparison was performed between these statistics before and after the stimulation. The results (Table II), show no significant difference in the average values of these statistics in the presence of TMS. The findings, therefore, suggest that in the studies related to the effectiveness of TMS in seizure-control, these parameters have no significant performance and if a PA is being designed to extract these features, its accuracy would not be acceptable. At the next step, the LLE was computed for every recording. The average values of LLE, are presented in Table III. The results revealed that in the ictal phase, the LLE was mostly larger than the pre-ictal phase. LLE had its least mean-value for the normal volunteers and its largest mean-value for epilepsy patients in their ictal phase. But, in 16 recordings, the LLE of the pre-ictal phase, was higher than its value in the ictal phase. Therefore, if the PA is defined just based on the LLE measures, its accuracy would be 73.3% which is not acceptable for clinical applications. Therefore, the detection of bifurcation was studied as the next step. For each recording, a specific threshold was selected with trial and error. Then the ISIs at different values were calculated for distinct windows of the data. The results are presented in Table IV. Bifurcation detection was carried out in the first, second, third, and forth windows. The accuracy of the detection was 100% in the pre-ictal phase and mostly between the first and second windows of data. The bifurcation diagrams for the recordings#9, 54, and 47 are illustrated in figures 6-8. For the pre-ictal recording#9, the bifurcation was detected between the first and second windows. Each window contains 1.15 seconds of 23.6 seconds of data. Hence, after the detection, there is at least approximately 21 seconds before the seizure onset, which should be considered in computing the final detection timelength. For the pre-ictal recording#54, the bifurcation was detected between the second and the third windows. The threshold voltage was calculated to equal 56 mV and the bifurcation was detected after 2.3 seconds. For the recording#47, the bifurcation was detected after 3.45 seconds and between the third and the forth windows. The average time necessary for computation of all the 60 non-real-time recordings was 0.3443 seconds.  The results of the real-time computation to design a realtime PA, are illustrated in Table V. The average training time was 2.19 seconds for all recordings. The training can be performed several days to several minutes before the application of the PA. To achieve a complete comparison between the real-time and non-realtime PAs, the sensitivity, specificity, and accuracy of the algorithms, should be calculated. In this study, the true positive detection (TPD) applies when the bifurcation is detected in the ISIs of pre-ictal recordings with respect to the induced magneto-motive force, . False-positive detection (FPD) applies when a bifurcation is detected in the ISIs of an ictal EEG recording. True negative detection (TND) applies when the bifurcation is not detected in the ISIs of the ictal phase. False-negative detection (FND) applies when the bifurcation is not detected in the ISIs of the pre-ictal phase. According to these definitions, sensitivity, specificity, and accuracy can be defined as [29]: Using these definitions, a comparison between the nonreal-time and real-time algorithms has been presented in table VI. According to the results, the sensitivity, specificity, and accuracy of the non-real-time PA, were better than the real-time one, but the prediction speed of the real-time algorithm was considerably higher than the non-real-time PA. Moreover, in the real-time algorithm, the training phase which is used to defining the ℎ was performed dramatically faster. The whole results suggest that the proposed PA is highly applicable in clinical applications.  [33] and one of which is operating realtime [34]. The features extracted in each algorithm are reflected in the table. Some of these researches such as [3] and [6] have performed the computations on the same dataset used in the present study. Comparative analysis of the sensitivity, specificity, and accuracy of the previous algorithms and the proposed algorithm in this study, reveals that both proposed PAs have the best performance.

IV. CONCLUSION
A real-time seizure prediction algorithm that is applicable to a non-invasive neuromodulator was proposed in this study. The main stone of the algorithm was detecting the bifurcation in the pre-ictal phase of the electroencephalograms of the epilepsy patients. A comparison between the sensitivity, accuracy, and specificity of the proposed algorithm and other prediction algorithms revealed a significant performance of the proposed algorithm in clinical applications.
One of the breaking points of the simulations in the current study was the lack of real EEG data under the TMS. To address this issue, a TMS-compatible EEG-recording device should be used to achieve the real EEGs under the TMS. Simultaneously, a Hall Effect sensor could be used to calculate the distribution of the magnetically induced voltages and currents.
To enhance the accuracy and sensitivity of the prediction algorithm, several actions can be addressed. Since it has been shown that the number of neurons of each functionally specified brain region is about 10 5 to 10 6 neurons, increasing the number of neurons of the SWN would enhance the similarity of the SWN to an actual epileptogenic zone.
Another task could be studying the effects of the current passing through the coil on the and . This would pave the way of using this current as a control parameter and design a responsive neuro-stimulator which can detect the bifurcation and select the appropriate coil's current to replace the starting point of the bifurcation and to extract a desirable response from the network. One of the controllers that can be proposed to suppress the seizure-like activity in an SWN, is the time-delayed feedback control system. By dividing the main system into linear and nonlinear parts, the periodic solution of the system and its frequency can be estimated. Regarding this frequency, the magnitude of time delay can be estimated too. To achieve an optimal control gain, the main system should be simulated and the feedback gain would be defined so that the controller could extract a stable periodic response among the chaotic responses.

A. Huber-Braun equations
Intended to study the effects of the magnetic stimulation on a small-world-network of Huber-Braun neurons, the first step should be simulating the coil of the stimulator. Since the circular coils can stimulate deeper structures [34], a circular coil with a radius of 1 cm and 15 wire turns were simulated.
Computing the , necessitated the knowledge of the induced electrical field, E, composed of two components: and . , is the time derivative of the magnetic potential, , and is related to the gradient of charge accommodation on the air-tissue interface, .
If the coil is located parallel to the air-tissue interface, then is negligible [35]: Simulation results revealed that the maximum value of ⃗ is just under the coil periphery. As the distance from the coil increases, there would be a less strong but more localized field. Based on the fact that the distance from the scalp to the cortex of an adult is approximately 2 cm [36], the optimal nerve-to-coil distance was considered equal to 2 cm. The HB model equations are as follows: = ( ) ( − ). : , , , ( ) is a temperature-dependent coefficient and is described as follows: ( ) is the activating coefficient, is the time constant and ∞ is the steady-state activating factor that is described by the following equation: is the steepness and 0 is the half activation value. The above equations are used to describe and , but is described as: is a complex current that is made mostly by Na ions and the K ions have less influence on it, but the latter ions are essential in producing the slow current. represents the calcium-dependent potassium current and can be computed when the dynamics of the internal changes are calculated by the equation of the activation dynamic, : where describes the increase in the concentration of +2 ions but is the amount of active cancellation of intercellular calcium [18]. Table I shows the values of the parameters of the aforementioned equations.

B. Huber-Braun Cable Model
The passive cable equation of a nerve fiber is described in Fig.1, where ( ) is the intracellular current over the fiber axis. ( ), is the intracellular potential and ( ), is the fiber's membrane current per its unit length.
is the membrane potential and , is the inserted voltage from outside of the cell. While the cable is inactive, the induced external voltage is equal to zero and the membrane voltage is equal to the internal potential: = , By combining these equations, the passive cable equation is achieved as: The above parameters are: = √ : the length constant = : the time constant, = 2 : the fiber resistance per unit length, : the resistance of the axoplasm, : the fiber's radius =2 : the membrane current per unit length, : the membrane's current density, =2 : the membrane capacity of unit length, : the membrane's capacity per unit area and : the membrane's resistance per unit length. Fig. 1. The electrical circuit of the passive cable model [35]

C. Modification of the Huber-Braun Cable Model
The active cable equation is obtained by applying TMS. It was assumed that and are the components of the electrical field, which are inside the nerve fiber and along its axis respectively: Research elsewhere show that has a more significant effect in depolarizing the membrane than two other direction components. It is on such a basis that their effects are negligible [36].
The modified active cable equation of a nerve fiber under the TMS is achieved by combining the above equations: The parameter, , is considered as the equivalent to the Rattay's activating function in magnetic stimulation [36].
Applying the induced elements to the classic HB model, the modified equation was obtained as follows: where 2 is the axial Conductance of the cable.
Since the parameter, , depolarizes the membrane, hence: The solution of this partial differential equation, necessitated the application methods of lines (MOL) to convert that to an ordinary differential equation. The methodology involved dividing a finite nerve fiber with the length, , into shorter length elements, ∆ , and the secondorder partial derivatives are converted to the difference elements [37], Therefore, (25) can be written as: According to the above formulation, external induced factors can be considered. This made it possible to investigate the different model behavior due to different measures of the induced current. On the other hand, the membrane potential is the most important variable to investigate different behavior of the classic HB model. The neural messages are transformed via the spikes of the membrane voltage. The inter-spike-intervals (ISIs) that are the distances between two consequent spikes are considered to be the most important variables for calculation [3][22] [23].

D. Small-World-Network of the Huber-Braun Model
After modification of a classic HB fiber, an SWN of the modified HB neurons was provided as a model of the epilepsy foci. The initial measures in a small-world-network are the average distance between its nodes, L, and the clustering coefficient, C [21]. The latter means the level of overlapping between two neighbors from two opposite sides. In other words, if the node , is connected to the nodes and , then C defines the probability that the nodes and are connected. The SWNs have small L and big C. Two parameters that describe the network's architecture are N and P. The former refers to the number of neurons in the network and the latter is the probability of nonlocal connections. P and N were assumed by [18] to be 0.001 and 2000 respectively. The normalized measures of C and L were, therefore, considered to be 0.9 and 0.2 respectively. These measures are thought to satisfy the requirements of the small-world-networks.
The aforementioned modified HB model was used to simulate the nodes of the network. The chemical synapses were modeled as a means of simulating the connections between the physically far apart neurons. The synaptic dynamics can describe the impulsive effects of a presynaptic neuron on a postsynaptic one, while the former is triggering an electric pulse [18]. The equation of the -th neuron of the network is: = ( ) ( − ), : , , , The inter-network connections are described by the adjacency matrix [ ], whose elements, , are 0, if neurons and are not connected and 1, if they are connected. , is the coupling strength and ( ), is the portion of the band-receptors of the -th neuron whose time evolution is described by the following differential equation: , is the membrane potential of the presynaptic neuron and and are the characteristic rise time and settling time of the chemical synapse.
The mean voltage of the SWN, , is a beneficial parameter to study network behavior.
When the neurons are not synchronized, represents a tonic-firing behavior. However, during the coupling phase, represents a bursting-spiking activity. To simulate the EZ, a small-world-network of 2000 modified HB neurons was provided. Considering as the bifurcation parameter, a bifurcation was observed in the ISIs of . The resulted diagrams are illustrated in the next section.

E. Simulation Results
The circular coil was considered as a 1000-sided polygon in the x-y plane as showed in Fig.2. With the assumption that the current amplitude being 3800 A and the frequency is 1 Hz, the distributions of the magnetic potential as well as the induced electrical field were computed and illustrated in Figures 3 and 4. Simulation results revealed that the induced electrical field has its maximum magnitudes at the distances near the coil's periphery and just under its projection on the x-y plane. Increasing the distance from the coil plane, resulting in a more focal and localized field, but its strength was being decreased.
For each coil geometry and specification, an optimal coilfiber distance could be defined by trial and error. In this study, the computations were conducted assuming that this optimal distance is equal to 2 cm. So, the distribution of the activation function along the fiber axis (x-axis) at this distance, was computed and illustrated in Fig.5. Previous studies have revealed that the straight nerve fibers are stimulated more easily at the peaks of the activation function [38]. Hence, the maximums of the activating function were detected and the initial point of the nerve fiber was located at the first optimal point. To find the distribution of the ionic currents and membrane potential, the nerve fiber of 400 µm length was divided into 40 segments of 10 µm. Then, the membrane potential and the induced current were calculated at different times and locations. Depending on the different measures of , the qualitative different behavior was observed in the location where the stimulation occurred. For instance, where = −0.68 2 or less, the model showed tonic-firing behavior ( Fig.6(a)). By increasing up to 1.385 2 , the voltage spikes showed a decrement, but the distance between each two consequent spikes showed an increment and the neuron represented a bursting-spiking activity ( Fig.6(b)). Increasing the resulted in subthreshold oscillations (Fig.6(c)) and steady-state ( Fig.6(d)). Around = 1.385 2 , the modified model represented chaotic behavior (Fig.7).