Analysis of the behavior of a seizure neural mass model using describing functions.

Neural mass models are computational nonlinear models that simulate the activity of a population of neurons as an average neuron, in such a way that different inhibitory post-synaptic potential and excitatory post-synaptic potential signals could be reproduced. These models have been developed either to simulate the recognized neural mechanisms or to predict some physiological facts that are not easy to realize naturally. The role of the excitatory and inhibitory activity variation in seizure genesis has been proved, but it is not evident how these activities influence appearance of seizure like signals. In this paper a population model is considered in which the physiological inter-relation of the pyramidal and inter-neurons of the hippocampus has been appropriately modeled. The average neurons of this model have been assumed to act as a linear filter followed by a nonlinear function. By changing the gain of excitatory and inhibitory activities that are modeled by the gain of the filters, seizure-like signals could be generated. In this paper through the analysis of this nonlinear model by means of the describing function concepts, it is theoretically shown that not only the gains of the excitatory and inhibitory activities, but also the time constants may play an efficient role in seizure genesis.

INTRODUCTION

Models are created in order to study complex and interesting phenomena. The main purpose of modeling and simulation are the analysis and the understanding of observed phenomena, testing hypotheses and theories, and predicting the system behaviors in space and/or time. Models are useful, because studying a model should not be more difficult than studying the phenomenon it describes. The physiological mechanisms that underlie the transition to seizure in epilepsy in humans are still largely unknown. Both experimental[123] and theoretical models[456789101112131415161718] have been developed in an attempt to mimic the conditions under which seizure activity occurs. Theoretical models are usually systems of equations or computer simulations.

Although it is clear that modulations in excitation and inhibition in synaptically connected cortical neuronal networks can lead to abnormal discharge patterns, from the physiological observations it is not known appropriately how these modulations occurs. But it is possible to show by investigating which model parameters are effective in producing characteristic activities of a seizure. Thus, a central requirement in theoretic approaches is to demonstrate the parameters whose variations can generate epileptiform activities. This may provide some hypothesis for neurophysiologists to test whether they are physiologically plausible or not. Hence, the task facing the model designer is analyzing the model behavior with respect to the parameters.

Brain activity recording and thus its modeling approaches have been considered at different levels: Cellular (microscopic), population of neurons (mesoscopic) to macroscopic. As known, the brain activity is the result of the neurons’ interactions. Therefore, the behavior of individual neurons and their interconnections have been considered at the microscopic level, in such a way that by setting suitable relations among thousands of neurons the local field potentials could be modeled. Modeling the individual neurons is the most precedent one in the field of expressing neuro-physiological phenomena by some mathematical and computational equations. Hodgkin-Huxley model[19] is one of the most famous models for the neuron. Noisy and leaky integrate-and-fire neurons and a Poisson spike-train cell model are the other examples. Specialized models for hippocampal neurons are also proposed.[1011] These models mainly differ in the nonlinear function of dendrite, the linear impulse response of the soma filter and modeling refractory phenomenon, but are common in being stochastic models. Netoff[20] modeled the brain as small world network in which all cells are only coupled to their nearest neighbors, but small numbers of connections are broken and rewired to make long-distance connections at random locations.

In neuronal tissue (gray matter of the brain) neurons can often be divided into two main families, the principal neurons and interneurons. Thus, the models of neurons have their special parameters for each kind of neuron according to their properties. On the other hand the neurons of a neural network can be connected to each other in different lattice forms. In these microscopic level models any change in the ratio of the number of excitatory to that of inhibitory neurons, as well as the structure of the connected neurons and the relevant coefficients can change the type of activity at the model output. In other words, there are too many effective parameters in these models that prevent the model to be suitably tractable. Thus generally, analysis of the microscopic models to discover some physiological facts is very difficult.

In another modeling scheme, called the reduced model, the complexity measure of single neurons is low so that a few neurons can build Electroencephalography (EEG) signals. In[21] cortex is modeled as a matrix that each component contains excitatory and inhibitory mass of neurons. These models are mainly modified versions of the microscopic models that can lead to better comprehension of cortex dynamics.

For gaining a better insight into the neural functions underlying the emergence of different dynamics, further simplification of the models have been proposed, in which neuronal networks of cortex have been considered as a spatially continuous network. In these macroscopic level models, the properties of a bulk of neurons at a local region have been described by time-space dependent state variables. Typically, these variables can be the mean firing rates or the mean value of the cell (soma) membrane potential, just like there is an average neuron at that location.[1822] These kinds of models are also called neural mass models. Analogues to a single neuron[6] nonlinear function at the average dendrite and linear function at the average soma is expected in a neural mass model. Furthermore, the input of a neural mass model is usually a white noise to account for the influence of all far neuron populations.

Although these models are the simplest ones introduced for cortex, their differential equations possess nonlinear terms, and thus are difficult to analyze. Almost no analytical solutions exist for these equations and they must be qualitatively analyzed. By assessment of the model[2324] further insight has been gained into the properties of the dynamics of cortex performance. The strength and time course of cortical inhibitory–inhibitory and inhibitory–excitatory connectivity have been realized to be the effective parameters of the model about the seizure initiation and cessation. Furthermore bifurcation diagrams have been considered to analyze the model behavior.[25] Also, by linearizing the part of the model responsible for the fast onset activity about its equilibrium solution, Molaee-Ardekani[26] quantified the characteristics of the resonance frequency of that part of model while free parameters vary. It has been shown that if the parameters are suitably selected, sustained limit cycle activity with high frequency and damped oscillatory behavior in response to brief stimulating inputs can be seen.[26]

Although linearization and the associated sensitivity analysis are able to provide a great deal of information regarding the conditions under which linear stability is lost, they are in general insufficient for a nonlinear model. In this paper the describing function approach is used to analyze a neural mass model more accurately, in order to show if other parameters of the model can influence the model output activity. To begin, we focus on a single parameter: The excitatory time constant of pyramidal neurons. It is concluded that as well as the excitatory and inhibitory gains, the excitatory time constant of pyramidal neurons has a significant role in initiation and termination of seizures.

MATERIALS AND METHODS

Depth-EEG Model

The first neural mass models only considered the excitatory interactions.[27] Later, Wilson and Cowan[2829] incorporated both excitatory and inhibitory properties of the neuronal populations connected to each other. The state variable was the mean firing rate, and the equations underlying these models are differential equations that can easily be extended into partial differential equations, to account for spatial properties of the neuronal mass. Amari[3031] has proposed a spatio-temporal model with the mean soma membrane potential as macroscopic state variable.

Lopes Da Silva et al.[6] and simultaneously Freeman[32] developed a lumped-parameter population model as a forward model, which is capable to explain the alpha rhythm of the EEG and deals with perceptual processing in the olfactory system. Physical laws they employed are related to two conversion operations involved in the dynamics of each neural ensemble: A mean wave-to-pulse operator at the soma of neurons, and a mean linear pulse-to-wave conversion implemented at a synaptic level.

Wendling et al.[14] used the same approach to interpret depth-EEG signal characteristics in epileptic patients according to physiological facts. Indeed, as a new law, they tried to model realistic epileptiform activity as a result of imbalance between excitatory and inhibitory synaptic gains (i.e., model parameters). In this model each area of hippocampus contains three interacting subsets of neurons. Main cells (i.e., pyramidal cells in the hippocampus or neocortex) receive excitatory and inhibitory feedbacks from other subsets composed of local interneurons. Also in order to represent fast EEG activity such as low-voltage rapid discharges that are often observed in depth-EEG signals at seizure onset, Wendling et al. developed a new version of mentioned model based on bibliographical material.[14] The second class of inhibitory interneurons is added to model as the fourth subset to raise kinetics of the system. In this local model the influence from neighboring areas of hippocampus is represented by an excitatory model input, x(t), which is represented by Gaussian white noise. Thus the model in[14] is a stochastic model whose output corresponds to the post-synaptic activity of the first subset, and is interpreted as depth-EEG signal. As mentioned before, in this paper we call it the depth-EEG model.

Liley developed the basic theories of the macroscopic models by adding some new physiological facts.[18] Thus, effects of the synaptic reversal potentials, which make the amplitude of the respective post-synaptic potentials dependent on the ongoing post-synaptic/somatic membrane potential, have been added to the model. The effect of synaptic reversal potentials is the changing of the linear relation between the mean soma membrane potential and synaptic input to an exponential one. Also transmission of axonal pulses in long range fibers has been added to this model as a function. In other words, Liley model is capable of reproducing the main features of spontaneous human EEG. In particular, autonomous limit cycle and chaotic oscillatory activity in the alpha band (8-13 Hz) has been easily produced.

As shown in Figure 1, in depth-EEG model each synaptic process (excitatory, fast and slow inhibitory processes) is a combination of a linear system (pulse-to-wave conversion implemented at a synaptic level) and a nonlinear function (wave-to-pulse operator at the soma of neurons, i.e., S(v) =2e0/(1+e)), where v is the average potential of the pre-synaptic cells and S (v) is the mean firing rate of post-synaptic cells. This process is expressed with the differential (1).

Figure 1

Organization of the depth-EEG model. Four neuronal subsets: pyramidal cells, excitatory interneurons (A), dendritic projecting interneurons with slow synaptic kinetics (GABAA,slow:B) and somatic-projecting interneurons (the grey rectangle) with faster synaptic kinetics (GABAA, fast:G). The average pulse density of afferent action potentials is changed into an average inhibitory or excitatory post-synaptic membrane potential using a linear dynamic transfer functions, while this potential is converted into an average pulse density of firing of post-synaptic neuron using a static nonlinear function (from [14])

The impulse response of the second-order linear systems involved in the model is characterized by the time-constant and the peak amplitude (that is reached after one time-constant), as shown in Figure 2. These are two indicators of the system. We call the former “synaptic time constant”, and the latter is proportional to the “synaptic gain”. It is evident that fast-inhibitory synaptic time is less than slow inhibitory synaptic time. According to[14] standard values of linear system (C1 to C2, a, b, g) and nonlinear function parameters (e0 v0 and r) are presented in Table 1[14] y0 − y9 are the states of the model, which are equivalent to the excitatory post-synaptic potential (EPSP) and inhibitory post synaptic potential (IPSP) of the subpopulations of the neurons, or their corresponding mean firing rates.

Figure 2

Average standard post-synaptic membrane potentials: excitatory, slow inhibitory and fast inhibitory, respectively, obtained from impulse responses given by he(t) = Aate−ta, hi(t) = Bbte−tb. and he(t) = Ggte−tg, t> 0 for A= 3.25, B= 22, G= 10 (from [14])

Table 1

Model parameters, interpretation and values used to produce EEG signals adapted[14]

Types of the Model Output Activity

Based on the assumption that at different vigilance and pathological states the physiological parameters of the population of hippocampal neurons may differ from standard ones, Wendling et al. have produced six different types of EEG signals by changing the excitation and inhibition synaptic gain parameters (A, B, and G, respectively) in the (0,100) interval, which is a realistic range.[14] Normal background, sustained discharge of spikes, low voltage rapid activity, slow quasi-sinusoidal activity, sporadic spikes, and slow rhythmic activity. Each type of these signals can be met in the real world (normal and pathological) during different activities of the brain.

Figure 3 shows typical output signals of model according to different parameters (i.e., different models of the model space). Since (1) are stochastic differential equations, they should be solved by a stochastic numerical algorithm; but, like[33] we accept the simplicity of solving our equations by Euler algorithm.

Figure 3

Six different kinds of the depth-EEG model output

We used Gaussian white noise as hippocampcal modulatory input. In order to produce background activity parameters using standard values,[14] the mean and variance of the input signal should be properly chosen. To this end, we selected μ = 50 and σ = 50.

Motivation to Analyzing Model

It has been observed that if all parameters of the depth-EEG model are supposed to have constant values equal to the standard ones except excitation/inhibition synaptic gain parameters (A, B, and G) some differences between the model output and real signals, especially in dominant frequency of low-voltage rapid activity, are obvious. It is seen that dominant frequency of low-voltage rapid activity in the model output is non-flexible, but in the real world it is not the case. In fact, frequency of the seizure activity is not unique for different subjects and also for different seizures of a person. Although these changes are small they should be considered.

To comply with the possible changes of frequency, some other parameters must be assumed to be variable over time. The question is then which parameters could be effective in this order. According to the simulation results synaptic time constant of pyramidal cell excitatory activity affects the frequency content of the model output. The block corresponding to pyramidal cell excitatory synaptic process of the model is bolded in Figure 1.

Thus analyzing the model is necessary to indicate how the time constant of pyramidal cell excitatory activity influences the model output frequency and even the activity type.

Describing Function Concept

A nonlinear system can be studied by simply replacing the nonlinear operations by approximating linear operations. When the operating point of the system is changing by time, instead of repeating the linearization about new operating points, one must determine the operation performed by the nonlinear element on an input signal of finite size and approximate this in some way by a linear operation. Thus the same nonlinearity when driven by inputs of different forms, or even when driven by inputs of the same form but of different magnitude, must be used. The approximation of a nonlinear operation by a linear one which depends on some properties of the input is called quasi-linearization.

Three basic signal forms with which the quasi-linear approximators for nonlinear operators have been driven until now are: Bias, Sinusoid, and Gaussian process. Also, different linear combinations of these basic input signals are considered. The quasi-linear functions that describe approximately the transfer characteristics of the nonlinearity are termed describing functions (DF).[3435] The concept is to see when, for example, a sinusoidal signal is applied to the nonlinear function, what would be the fundamental component of the output signal. The ratio of the (output/input) phasors at the fundamental harmonic is defined as DF of the nonlinear function for sinusoidal input. This concept is mainly used in the detection of the system's limit cycle, i.e., prediction of its approximate amplitude and frequency, to be able to design suitable controllers for nonlinear systems.

To obtain the DF for a sinusoidal input, x[n] = a0 sin[ω0n], it must be considered that at steady state the output of the nonlinear characteristic, y[n], is periodic but, in general non-sinusoidal. Assuming the nonlinearity to be symmetric about zero, the Fourier series becomes:

where Y and φ are the amplitude and the phase shift of the kth harmonic component of the output y[n], respectively. In the sinusoidal-input DF analysis, only the fundamental harmonic component of y[n], Y1 is considered, because the higher-harmonics are often of smaller amplitude than the amplitude of the fundamental component.

Accordingly, the sinusoidal-input DF of a nonlinear element, N(a0,ω0), is defined as the output/input ratio of the fundamental harmonic component.

DESCRIBING FUNCTION OF THE MODEL's NONLINEAR BLOCK

To understand accurate role of the depth-EEG model variables the nonlinear blocks can be approximated by a linear function equal to their suitable DF.[3435] There are four nonlinear blocks, which in Figure 1 are labeled by numbers 1-4: N1, N2, N3, N4.

When the periodic activity is seen at the model output, input of the nonlinear blocks can be assumed to be sinusoidal signals with a bias, such that “sine-plus-bias DF” would be sufficient to analyze the model. Since the nonlinear function S(v) does not cause any delay, DF has no phase-shift term. So only the magnitude of the output's fundamental harmonic and the bias value must be determined. Accordingly, to obtain the DF for S(v) the integrations (2) must be computed:

where a0, ω0 and y0 are amplitude, frequency and the bias value of input signal of the nonlinear block. Moreover, a0 and y0 are the amplitude and bias of the output signal of the nonlinear block.

As mentioned before, the nonlinear function of depth-EEG model is S(v)=2e0/(1+er(v). So, two integrations of (4) would be:

In order to simplify the above integrations to be analytically calculated, the nonlinear function is approximated by a piecewise linear function. Since the sigmoid function S(v) is a smooth and strictly monotonic function, the piecewise linear approximation is a convenient one. Although the larger the number of these linear pieces, the better the nonlinear function would be reconstructed; but, using only three linear pieces suffice in our analytical calculations. This is firstly because the approximated nonlinear functions do not lead to alter the model outputs significantly, and secondly the effect of other parameters of the model, which is the goal of our analysis, will be evident without any necessity to take the complexities into account. The approximated function is shown in Figure 4, and its related equation is displayed in (6).

Figure 4

The non-linear function is approximated by a piece-wise linear function

Accordingly, the 0 – 2π interval of integration must be partitioned into smaller intervals, in which S (v) is linear. For example, to this end in Figure 5 θ1 and θ2 indicates the borders of sub-integrations. Figure 5 shows the input and output of the nonlinear function S (v) where the sinusoidal amplitude equals to a0 = 16 and its bias value is y0 = –2.5.

Figure 5

Example of the time intervals used in order to integration, according to the approximated nonlinear function

It is evident that depending on the values of a0 and y0 different situations occur, for which the integration intervals must be indicated separately

Before describing these different situations, note that θ1 and θ2 are equal to and respectively. Since x1 < x1 then θ1 < θ2. Furthermore, due to the symmetry of the sinus signal, both θ1 and θ2 have values between –π/2 to π/2. Thus, three different cases have been detected.

Case 1: θ1 < 0,θ2 < 0, or equivalently , thus

Special cases of this situation are:

Case 1a: –π/2=θ1, θ2 < 0 that occurs when and :

Case 1b: –π/2=θ1,–π/2=θ2, that occurs when , then:

Ŝ(y0+a0sin[ω0n+α]) = 2e0      (9)

Case 2: θ1 < 0, θ2 > 0, or equivalently and

that its special cases are:

Case 2A: –π/2=θ1,θ2 > 0, or equivalently and

Case 2b: θ1 <0,θ2 = π/2, or and

Case 2c: θ1 = –π/2, θ2 = π/2, or and

Case 3: 0 < θ1, 0 < θ2, or equivalently

with the special cases:

Case 3a: θ1 >0, θ2 = π/2, or

Case 3b: θ1 = π/2, θ2 = π/2, or

Ŝ(y0+a0sin[ω0n+α]) = 0      (16)

After computing the integration and some algebraic operations, the result of all of these cases for the values of 0 and 0 can be written in the same closed form as following:

For example for case 2c that θ1 = –π/ 2,θ2= π/ 2, the output of the S(v) in response to a sinusoidal input with amplitude a0 and the bias value y0, is a periodic signal whose fundamental harmonic amplitude is 0 and the dc component have the following values:

ANALYSIS OF DEPTH-EEG MODEL AT A SPECIAL PARAMETER VECTOR A = 21, B = 31 and G = 21

To investigate whether other variables of the depth-EEG model are effective in changing frequency contents of the model output and/or initiation or termination of the periodic seizure-related activity, values of A, B, G have been fixed. The problem is limited to assess the role of a single parameter, here of the pyramidal excitatory time constants (τ). We want to show the validity of the following:

Periodic activity can appear or disappear by changing the value of the pyramidal excitatory time constants (τ).

Frequency of the periodic activity is dependent on the τ.

Although detail of the analysis approach is dependent to the values of A, B, G, but the main approach is the same for all parameters: Different assumptions about the cases of nonlinear blocks should be tested until valid cases for all range of τ would be found. Thus, in this paper to describe how the analysis would be done in details, as an instance, a special value of A, B, G is considered: A = 21 B = 31 G = 21. These values are manually selected and there is a periodic activity for a = 100, b = 50, g = 500, τ = 50 at the model output. To obtain a general result, the depth-EEG model could be analyzed for other values of A, B, G. Also, the role of other time constants could be indicated.

The Linear Part of the Depth-EEG Model

The linear filters of the model are low-pass filters. In simulation with MATLAB software they are treated as discrete filters by the following transfer functions.

where T is the sampling time by which the differential equations are numerically solved. The dc gain of filters (ω = 0) will be A/τ, B/b, G/g, and A/a. The frequency responses of these filters are as follows:

The Nonlinear Part of the Depth-EEG Model

According to the monotonic nature of nonlinear functions, it is clear that at every set of the model parameters there are unique state signals and output signal for a special input signal. Thus, depending to the dc value and the amplitude of the fundamental harmonic of the nonlinear blocks, they may work just in one of the abovementioned cases. To determine the valid cases of under the assessment parameters, we start our analysis by assuming some test fixtures for the nonlinear blocks. By analysis of the equivalent linear system, the dc value and the amplitude of the fundamental harmonic of the model output signal can be easily computed. By using these values, we obtain the input of every nonlinear block and check it to see whether or not the assumed case for that block is correct. If the assumed cases were not the correct ones, i.e., leads to some contradictions in the system, we must continue the analysis procedure by assuming some new cases. Once the valid set of cases for the model parameter is found, the analysis is finalized for that parameter.

We have checked different assumptions by trial and error, but here just those assumptions that have been found to be correct at least for some values of τ in under the assessment range (0,100), are reported. Thus, at first we suppose that for A = 21, B = 31, G = 21, τ = 50 , a = 100, b = 50, g = 500, N1 is at case 2c, N2 and N3 are at case 1b, and N4 is also at case 2c.

Assumption 1: N1 and N4 are at Case 2c, N2 and N3 are at Case 1b

Assuming the nonlinear blocks of the model at the above mentioned cases, the output of N2 and N3 are equal to 2e0. In other words, according to Figure 6, y1 and y2 are, respectively, obtained by passing a dc value (2e0) multiplied by C2 and C4, through filters H (e) and H (e). Thus, y2 and y2 are constant signals, respectively, equal to 2C2 Ae0/a + μ and 2C4 Be0/b, where μ is the mean value of model input noise. Under these assumptions the effective loop of the model in which the oscillation is expected to be generated is shown in Figure 6.

Figure 6

The approximated model by using describing function under Assumption

Furthermore, since N1 is supposed to work in Case 2c it would be replaced by its describing function, i.e., a gain equal to re0/In (99) (18). Similarly N4 is replaced by the same gain.

Assuming the y0 to be the notation for the dc value of the model output (input of the N1), and a0 for its amplitude at the fundamental frequency ((y0+a0sin[ω0n+α])), also by defining and as the outputs of N1 and N4 respectively, the (18) can be used to obtain:

and

Accordingly, the dc loop equation

can be written as:

Since we are analyzing the effect of changing the value of τ on the system in an special case in which A = 21, B = 31, G = 21, a = 100, b = 50, g = 500, by substituting these values in the (23) the amount of y0 versus τ is obtained as follows:

On the other hand using the Mason rule for the loop, it is evident that the model is a simple linear system with the frequency response as follows:

Because the model input is a Gaussian white noise, the power spectrum of the model output by the assumption 1 about the situations in which the nonlinear blocks work would be as follows:

Checking Assumption 1

It is noteworthy that for some values of τ the (24) may not hold. In other words it must be checked when the conditions for Assumption 1 hold. The variance of the noise is assumed to be negligible in such a way that dc values of the state signals of the model are sufficient to check the conditions.

To check whether or not the above assumption about cases of nonlinear blocks is true, the following conditions must be assessed.

To check if N1 is in case 2c, it must be checked whether

To check if N2 is in case 1b:

To check if N3 is in case 1b:

To check if N4 is in case 2c, it must be checked whether: .

The above statements would be reduced to the following ones by some simple algebraic manipulations:

N1 is in case 2c if: 1.9774 < τ < 90.9735.

N2 is in case 1b if: 9.1190 < τ.

N3 is in case 1b if: 2.9605 < τ.

N4 is in case 2c for all values: 0 < τ < 100.

Totally the condition for which the above analysis is valid is as following:

9.1190 < τ < 90.9735      (27)

In other words the (26) is just valid for those values of τ with the limits of (27). This filter is equivalent to a resonant low-pass filter, which according to it a periodic activity with the resonant frequency is expected at the model output. The magnitude of frequency response (26) is plotted in Figure 7 for some values of τ in the interval of (27).

Figure 7

Magnitude of the frequency response of (26) for some values of t in the interval: 9.1190 < τ < 90.9735

Assumption 2: N1 is at Case 2c, N2, N3 and N4 are at Case 1b

Since in practice we consider τ to be in the (0, 100] interval, the system must be analyzed again for 0 @ τ < 9.1190 and 90.9735 < τ < 100. In this section the analysis is done for 0 < τ < 9.1190 by another set of assumptions about the cases of the nonlinear blocks. We now suppose that N4 also works in the case 1b, i.e., its output is equal to 2e0.

For Assumption 2, there is no loop in the system such that for dc we have the following equation:

which leads to y0 = 0.7665 for μ = 50. In this situation the model output is just a filtered white noise whose power spectrum is equal to:

We must check again the range of τ for which this situation occurs. It is noteworthy again that the variance of the input noise is supposed to be sufficiently small such that just the dc values indicate the cases for each nonlinear block.

To check if N1 is in case 2c, it must be checked whether

To check N2 is in case 1b:

To check N3 is in case 1b:

To check N4 is in case 1b:

By the obtained value y0 = 0.7665, (1) always holds, (2) requires τ ≤ 98.5814, (3) requires τ ≤ 24.6459, and (4) requires τ ≤ 9.1239, which all of them hold under the assessment range 0 < τ < 9.1190.

The spectrum of (29) is clearly a first order low-pass filter such that no periodic activity is seen at the model output for 0 < τ < 9.1190.

Assumption 3: N1, N2 and N3 are at Case 1b, N4 is at Case 2c

For 90.9735 < τ ≤ 100 the analysis is continued by assuming the case 1b for N1, N2 and N3, and case 2c for N4. In this situation too, there is no closed loop, and the equation for dc value is as the following:

that leads to:

For this situation, the spectrum of the model output is as follows:

Again by assuming the small variance for the input noise, the range of suitable τ for this situation is obtained.

To check if N1 is in case 1b:

To check if N2 is in case 1b:

To check if N3 is in case 1b:

To check if N4 is in case 2c, it must be checked whether:

These conditions lead to 94.7794 < τ < 119.1861. In other words for 94.7794 < τ < 100 the (32) presents the valid power spectrum of the model output. This means that for 94.7794 < τ < 100 no periodic activity is expected at the model output. For 90.9735 < τ < 94.7794, the analysis must be done again.

Assumption 4: N1 is at Case 2a, N2, N3 and N4 are at Case 1b

For 90.9735 < τ < 94.7794 we assume that N1 is at case 2a, N2, N3 and N4 are at case 1b. N1 behaves like a linear function by the following gain (θ2 > 0):

The spectrum of the model output would be very similar to 26:

This filter is again a resonant low-pass filter.

It is noteworthy that while A = 21, B = 31, G = 21, a = 100, b = 50, g = 500, by the mentioned assumptions the whole interval of [0,100] is covered. In other words, it is an easy task to show that any other possible assumption about the cases leads to a contradiction (incorrect result).

RESULTS

Based on the results of previous sections, it is evident that according to the concept of describing function, the complex function of the nonlinear model is approximated by a filter at A = 21, B = 31, G = 21. This filter is low-pass for 94.7794 < τ < 100 and < τ < 9.1190, but a resonant low-pass filter for 9.1190 < τ < 94.7794. In other words the role of the τ in changing the system function is now evident. But from the seizure initiation point of view, it must be discussed at which times the periodic activity appears in the model output. In Figure 8 changing the frequency response of the model by τ is shown.

Figure 8

Magnitude of the frequency response of the approximated model for A = 21, B = 31, and G = 31 for some values of t low-pass filter, the quality factor of the filter also must be considered to indicate whether the model output has a periodic activity or not. The quality factor of the resonant filter versus the value of is plotted in Figure 9 for 9.1190 < t < 94.7794

For the low-pass filters, 94.7794 < τ < 100 and 0 < τ < 9.91190, clearly there is no periodic activity in the model output, i.e., the normal activity. But for the resonant low-pass filter, the quality factor of the filter also must be considered to indicate whether the model output has a periodic activity or not. The quality factor of the resonant filter versus the value of τ is plotted in Figure 9 for 9.1190 < τ < 94.7794.

Figure 9

The quality factor values of the approximated resonant filter for the depth-EEG model, versus t, for. 9.1190 < t < 94.7794. In this figure the 0 < t < 9.1190 and 94.7794 < t< 100 regions for which the model is equivalent to low-pass filter are specified. According to the quality factor values for 90.3375 < t < 94.7794 is such small that the periodic activity of the model output is diminished

It is evident from Figure 9 that the quality factor of the resonant filter is decreasing as τ is increasing. The quality factor of the filter is so low for 90.3375 < τ < 94.7794 such that periodic activity seems to be diminished.

As a result we can say that at the constant value of all parameters, for example for A = 21, B = 31 G = 21, a = 100, b = 50, and g = 500, variation of the pyramidal excitatory constant time can initiate or/and terminate the periodic activity (seizure) by itself.

Also, the resonant frequency of the filter versus τ is shown in Figure 10. According to this figure variation of the dominant frequency of the model output as τ changes is obvious. Thus the pyramidal excitatory time constant is able to modify the frequency content of the model output. As it was mentioned before, unlike real signals, depth-EEG model suffered from its unique output seizure frequency that could be eliminated by taking τ as a new model parameter.

Figure 10

The resonant frequency of the approximated resonant filter for the depth-EEG model, versus τ, for 9.1190 < τ < 94.7794. In this figure the 0 < τ < 9.1190. and 94.7794 < τ< 100 regions for which the model is equivalent to low-pass filter are specified

These theoretical results are matched with that of simulation results.

DISCUSSION AND CONCLUSION

In this paper we showed the effects of changing the parameter τ for one special value of A, B, and G (A = 21, B = 31, and G = 21). However, the derived results can be generalized for other values. We can say that the excitatory time constant of the pyramidal neurons is another parameter that affects the frequency of seizure-related signals, i.e., decreasing the value of τ reduces the output frequency, and is also able to initiate and terminate seizure. That is, while a periodic activity is seen on the model output, it may disappear by changing τ (depending on synaptic gains), and vice versa: Sometimes changing τ may cause the appearance of rhythmic activity.

In other words, in the absence of any variation of synaptic gains, pyramidal cell excitatory time (τ) is able to induce or stop instability. This means that without any need for synaptic gains to return to their standard values, seizures could be stopped by changing parameter τ. This observation is in accordance with what Liley reported on his attempt[22] to find mechanism of both formation and cessation of generalized seizure activity, and propagation of focal epileptiform activity. As Liley declared, changes in fast GABAergic neurotransmission are necessary to induce transition to epileptiform activity, but it can be done through different parametric context. Weakening inhibition can be encountered by both reducing peak amplitude of the IPSP and its decay time constant. In this paper we showed that this may be the fact also for the excitatory processes: Every change in excitation process, either through change of the gain (A) or excitatory time constant τ can cause the instability or stop it. It proposes that the effective time of the excitatory and inhibitory processes, not only their gains, are important in changing the type of EEG activity. In fact integrating the excitatory and inhibitory processes seems to be more meaningful when they are compared with investigate if they are balanced or not. This fact must be assessed from the point of view of physiological mechanisms. The physiological meaning of changing the time constant of pyramidal neurons may resemble the change of the speed of neurotransmitters’ secretion at the synapses.

In this paper we focused on the effects of the excitatory activity of pyramidal neurons, but it deserves to evaluate the effects of other time constants in the model, e.g. a, b, and g.